Fast low-rank solution of the Poisson equation with application to the Stokes problem

E.A. Muravleva I.V. Oseledets

(August 8, 2025)

1 Introduction

Solvers for the Poisson equation are important in many different application areas, thus their improvement is of great interest even for special grids and right-hand sides. In this paper we consider a very specific class of Poisson equations: Poisson equations in two and three dimensions on tensor-product uniform grids with low-rank right-hand sides. In this case, it is possible to reduce the complexity of the solver to an almost linear complexity in a one-dimensional grid size. The problems of such kind have been considered previously by several authors [9, 6, 3] and the linear complexity is indeed possible. However, in small dimensions (especially in two dimensions) the constant hidden in $\mathcal{O}(n)$ can be very high, and the full representation is more efficient for a wide range of $n$ . The main goal of this paper is, staying within the framework of low-rank approximations, provide a new and efficient algorithm for the approximate solution the Poisson equation.

To show the effectiveness of our approach, we use it to create a fast solver for Stokes problem. The Stokes problem is one of the classical problems of mathematical physics. It is often encountered as a subproblem in more complicated problems, such as Navier-Stokes problem, unsteady Stokes problem, flows of non-Newtonian fluids. Numerical methods for the solution of the Stokes problem are a classical topic of the computational fluid dynamics, thus any improvement of the efficiency of the Stokes solvers (at least in some important cases) is of great interest. The Uzawa method for the Stokes problem requires the repetitive solution of the Poisson equations. The whole iterative procedure (including matrix-by-vector products, arithmetics, dot products) will be implemented in the low-rank format, and it is not always an easy task, since at each step the accuracy should be monitored to avoid the growth of the number of parameters and the growth of the error.

The method proposed in this paper has its limitations: it works for a special discretization of the Laplace operator on tensor-product grids and for special right-hand sides, but the class of right-hand sides in two and three dimensions is not small: it includes sufficiently smooth functions (i.e., approximated by polynomials), so-called asymptotically smooth functions, sparse right-hand sides and so on. The approximation work is done on the algebraic level using unified algorithms. To get better complexity and accuracy, one can use advanced discretization techniques: $hp$ -methods, high-order schemes, spectral and pseudospectral approaches, etc. However, they require significant additional work on the analytic level. On the other hand, one can use certain approximations on the discrete level, by approximating the discrete solution using some low-parametric format. As such a format, we will use a low-rank factorization of matrices (in two dimensions) and tensors (in three dimensions). The approach can be used in arbitrary dimension, but we will leave this topic for future work.

2 Model problem: discretization and solution method

The Stokes problem is used a basic problem to illustrate the low-rank techniques that are used. The numerical scheme consists of three basic steps.

1.

Take a simple discretization of the Stokes problem (we will use semi-staggered grids)
2.

Use a mesh-size-independent convergent iterative scheme (we will use the Uzawa method)
3.

Replace each step of the method by operations in the low-parametric format with truncation.

The whole procedure is very simple, and as it will be seen later on, most of the steps are also simple in the tensor format. However, to get a more efficient method, several modifications should be made to the “naive” approach. For the Stokes problem, the second step requires the solution of the Poisson problem at each iteration step. Our main result is a new algorithm for the solution of such problem, based on the cross approximation [23, 1] in the frequency space. The well-known approaches for the solution of the Poisson equation in low-rank formats rely on a special approximation to the inverse of the Laplace operator by a sum of tensor (Kronecker) products, but in two and three dimensions the complexity can be reduced signicantly using a new approach.

We consider the Stokes problem in a rectangular domain $\Omega=[0,1]^{d}$ , $d=2,3$ . The problem is discretized using the finite difference method (FDM) on semi-staggered grids. For details, see [16, 15, 14]. We will give only the final matrix form of the discrete problem. There are $d+1$ unknown vectors ( $d$ for velocity components and one for pressure). The components of velocity are defined in the vertices of the grid cells, and the pressure is defined in the centers of the cells (see Figure 1). The grid for pressure is $n_{x}\times n_{y}$ in two-dimensional case, and $n_{x}\times n_{y}\times n_{z}$ in three-dimensional case, respectively. For simplicity we will assume that the grid sizes are equal: $n_{x}=n_{y}=n_{z}=n$ . The mesh size is defined as $h=\frac{1}{n}$ .

Figure 1: The model problem: semi-staggered grids for velocity components and pressure

The discretized Stokes problem has the following form:

\begin{pmatrix}A&B\\ B^{T}&0\end{pmatrix}\begin{pmatrix}u\\ p\end{pmatrix}=\begin{pmatrix}f\\ g\end{pmatrix}.

(1)

We use the standard notations: $u$ is the velocity vector, $p$ is the pressure. Here $A$ is a $d\times d$ block diagonal matrix of the form

A=I_{d}\otimes\Delta,

where $\Delta$ is a discrete Laplace operator which will be defined later on, and $I_{d}$ is an identity matrix of order $d$ .

The matrix $B$ is a block matrix matrix of the form

B=\begin{pmatrix}B_{x}^{\top}&B_{y}^{\top}\end{pmatrix}^{\top}

in two-dimensional case, and

B=\begin{pmatrix}B_{x}^{\top}&B_{y}^{\top}&B_{z}^{\top}\end{pmatrix}^{\top}

in three-dimensional case.

To define the components of the discrete gradient operators, it is convenient to use the following auxiliary matrices, $G$ and $H$ .

They are defined as

G=E-Z,\quad H=E+Z,

where $E,Z$ are $(n-1)\times n$ matrices. The matrix $E$ deletes the first element of the vector (and moves all others up by one), and the matrix $Z$ deletes the last element of the vector (and moves all others down by one).

In two-dimensional case, the matrices $B_{x}$ and $B_{y}$ have the form

B_{x}=\frac{1}{4h}G\otimes H,\quad B_{y}=\frac{1}{4h}H\otimes G.

In three-dimensional case, the matrices $B_{x}$ , $B_{y}$ and $B_{z}$ have the form:

B_{x}=\frac{1}{4h}G\otimes H\otimes H,\quad B_{y}=\frac{1}{4h}H\otimes G\otimes H,\quad B_{z}=\frac{1}{4h}H\otimes H\otimes G.

Finally, the discrete Laplace operator has the following form:

\begin{split}&\Delta=B_{x}B^{\top}_{x}+B_{y}B^{\top}_{y},\\ &\Delta=B_{x}B^{\top}_{x}+B_{y}B^{\top}_{y}+B_{z}B^{\top}_{z},\end{split}

in two and three dimensions respectively.

Let us note that this is a non-standard 9-point stencil in 2D (27-point stencil in 3D) discretization for the Laplace operator, but its main benefit is that it is consistent with the discrete gradient operator.

The standard approach to solve (1) is the Uzawa method [2]. The system is reduces to the following equation for the pressure:

(B^{\top}A^{-1}B)p=B^{\top}A^{-1}f.

(2)

The matrix of the system (2) is positive semi-definite, and the conjugate-gradient (CG) method is a convenient tool to solve it. An important question is the spectrum of the Schur operator. It is known [16, 15] that it has $2$ zero eigenvalues in 2D, and $3n-1$ eigenvalues in $3d$ , and $(n-2)^{d}$ eigenvalues, equal to $1$ . The numerical experiments confirm, that the remaining part of the spectrum lies on $[0,1]$ , and is bounded from below by a constant, independent of $h$ . The linear system (2) is consistent, i.e., the right-hand side is orthogonal to any vector $q$ in the kernel of the Schur complement, thus the conjugate gradient method effectively works on the subspace, orthogonal to the kernel, where the operator is well-conditioned, thus the total number of iterations does not depend on $h$ .

3 Low-rank formats

Our main goal is the solution of the problem (2) in low-rank formats. In two dimensions, the discrete pressure is a vector of length $(n-1)^{2}$ , and it can be naturally considered as an $(n-1)\times(n-1)$ matrix. In a three-dimensional case, the discrete pressure can be naturally considered as an $(n-1)\times(n-1)\times(n-1)$ three-dimensional tensor. We want to approximate those tensors, effectively reducing the number of parameters and the computational cost, while maintaining the required accuracy $\varepsilon$ of the computations. This accuracy should be consistent with the mesh discretization error, which is $\mathcal{O}(h^{2})$ for the pressure.

Let us recall some basic facts about low-rank tensor decompositions in two and three dimensions.

Definition 3.1.

The $n_{1}\times n_{2}$ matrix $M$ is said to be in the low-rank format with rank $r$ , if it can be represented as

M=UV^{\top}=\sum_{\alpha=1}^{r}u_{\alpha}v_{\alpha}^{\top},

(3)

where $U$ is an $n_{1}\times r$ matrix, and $V$ is an $n_{2}\times r$ matrix.

The dyadic decomposition (3) can be computed via the singular value decomposition (SVD). The linear operator, acting on a space of $n_{1}\times n_{2}$ matrices, can be represented as a $(n_{1}n_{2})\times(n_{1}n_{2})$ matrix. The corresponding low-rank format for the matrix is given in the following form.

Definition 3.2.

An $(n_{1}n_{2})\times(n_{1}n_{2})$ matrix $A$ is said to be in the low-rank (Kronecker) format with rank $R$ , if it can be represented as a sum of Kronecker products:

A=\sum_{\alpha=1}^{R}C_{\alpha}\otimes D_{\alpha},

(4)

where $C_{\alpha}$ is an $n_{1}\times n_{1}$ matrix, and $D_{\alpha}$ is an $n_{2}\times n_{2}$ matrix, and $\otimes$ is a Kronecker product of matrices.

The low-rank format for operators and vectors allows for fast linear algebra operations [18, 19, 11]. For example, multiplication of a matrix of tensor rank $R$ by a vector of rank $r$ yields a vector of rank $Rr$ .

An important operation is the so-called rounding. Most basic operations (like addition or matrix-by-vector product) increase the rank, thus the reapproximation with some specified accuracy is required. If the matrix is given in the low-rank format, then the rounding procedure can be implemented without the computation of the full SVD. Its complexity with respect to $n,m,r$ is $\mathcal{O}((n+m)r^{2}+r^{3})$ , i.e., it is linear in the mode size (compared to cubic for the full SVD).

In three dimensions, the so-called Tucker format can be very useful: [22, 4, 13].

Definition 3.3.

An $n_{1}\times n_{2}\times n_{3}$ three-dimensional tensor $A$ is said to be in the Tucker format, if it can be represented as

A(i_{1},i_{2},i_{3})=\sum_{\alpha_{1},\alpha_{2},\alpha_{3}}G(\alpha_{1},\alpha_{2},\alpha_{3})U_{1}(i_{1},\alpha_{1})U_{2}(i_{2},\alpha_{2})U_{3}(i_{3},\alpha_{3}),

(5)

where the numbers $\alpha_{k}$ vary from $1$ to $r_{k}$ .

The matrices $U_{k}=[U_{k}(i_{k},\alpha_{k})]$ are called factors of the Tucker decomposition, the numbers $r_{k}$ are called Tucker ranks, and the three-dimensional tensor $G$ is called the Tucker core.

The Tucker format for the matrix is defined in the same fashion. The basic linear operations for the Tucker format can also be defined, as well as a fast rounding procedure with complexity that is linear in the mode size [18].

4 Fast inversion of the Laplacian in tensor formats

4.1 Known approaches

The basic computational core of the algorithm is the solution of the Poisson equation (see (2)):

\Delta u=g.

(6)

We will consider first the two-dimensional case. The generalization to higher dimensions will be outlined further. The matrix $\Delta$ in our case has the form (up to a scaling factor, which is not important)

\Delta=A_{1}\otimes A_{2}+A_{2}\otimes A_{1},

where

A_{1}=GG^{\top},A_{2}=HH^{\top}.

It is not difficult to see, that the matrices $A_{1}$ and $A_{2}$ commute and are diagonalized by the discrete sine transform (DST):

A_{1}=S\Lambda_{1}S^{\top},A_{2}=S\Lambda_{2}S^{\top},

thus the full matrix $\Delta$ can be written as

\Delta=(S\otimes S)(\Lambda_{1}\otimes\Lambda_{2}+\Lambda_{2}\otimes\Lambda_{1})(S^{\top}\otimes S^{\top}).

(7)

The representation (7) is a textbook way to solve such kind of problems in the full format, by applying the DST, inverting the diagonal matrix, and then applying the DST back.

For the low-rank formats, the situation is as follows. The application of the DST can be done in a fast way by computing the DST of the factors $U$ and $V$ , and it is an $\mathcal{O}(nr\log n)$ complexity. Moreover, then the rank does not change. But the inversion of the diagonal matrix

\Lambda=\Lambda_{1}\otimes\Lambda_{2}+\Lambda_{2}\otimes\Lambda_{1},

(8)

suddenly becomes the bottleneck, since it is a $\mathcal{O}(n^{2})$ computation! Thus, the matrix $\Lambda^{-1}\widehat{g}$ should be approximated in the low-rank format.

A well-established approach to compute the inverse of $\Lambda^{-1}$ is based on the approximation by exponential sums [21, 6, 9, 12, 10]. It is uses the following quadrature formula:

\frac{1}{x}=\int_{0}^{\infty}e^{-px}dp\approx\sum_{\alpha=1}^{R}w_{k}e^{-p_{k}x},

(9)

where $p_{k}$ and $w_{k}$ are quadrature weights. Typical value of $r$ required to achieve a good accuracy is of order of several tens.

The equation (9) can be used to get the approximation of the inverse matrix $\Lambda$ :

\Lambda^{-1}\approx\sum_{\alpha=1}^{R}w_{k}e^{-p_{k}\Lambda}=\sum_{\alpha=1}^{R}w_{k}e^{-p_{k}(\Lambda_{1}+\Lambda_{2})}\otimes e^{-p_{k}(\Lambda_{1}+\Lambda_{2})}.

(10)

4.2 A new method based on cross approximation

Let us estimate the typical complexity of the method, discussed in the previous subsection. The main computational task is the multiplication by $\Lambda^{-1}$ . Let $r$ be the rank of the vector, and $R$ be the tensor rank of the inverse of $\Lambda$ . Then the result would have the rank $Rr$ . Thus, in two dimensions, the method will be effective only when $n\geq(Rr)\sim 1000$ , since the typical values of $R$ and $r$ are around $30$ .

The main cost is due to the approximation of the full inverse matrix. However, we need only to multiply this inverse by a vector, and the result is expected to have a small rank. This structure is not used in the methods based on the approximation of the inverse matrix.

To obtain a new method, let us look more closely on the main operation:

\widehat{f}=\Lambda^{-1}\widehat{g},

or elementwise:

\widehat{f}_{ij}=\frac{\widehat{g}_{ij}}{\Lambda_{ij}}.

In our case,

\Lambda_{ij}=(4-\lambda_{i})\lambda_{i}+(4-\lambda_{j})\lambda_{j}=\mu_{i}+\mu_{j},

therefore

\widehat{f}_{ij}=\frac{\widetilde{g}_{ij}}{\mu_{i}+\mu_{j}}.

(11)

The most important observation is that (11) is an elementwise division of a matrix of rank $r$ by a matrix of rank $2$ , and we can compute any element of the matrix $\widehat{F}=[f_{ij}]$ in $\mathcal{O}(r)$ operations.

Our main assumption is that the result of the division can be in turn approximated by a matrix of a small rank, $\mathcal{O}(r)$ . Thus, we can try to approximate the result directly, without the approximation of the $\Lambda^{-1}$ . Since we can compute any prescribed element fastly, it is natural to use the cross method for the approximation of low-rank matrices (in two dimensions) [23, 1] and for tensors [17, 20]. In the matrix case, to approximate a matrix with rank $r$ it is sufficient to compute $r$ rows and $r$ columns of it, using a certain heuristical technique (a quasioptimal choice can be based on the maximal-volume submatrix [8] using the maxvol-algorithm [7]). Thus, the expected complexity is $\mathcal{O}(nr^{2})$ operations, which is much less, than the approximation of the full inverse.

5 Stokes solver in the low-rank format

The Uzawa method requires the solution of the Poisson equation and multiplication by the gradient matrices, and some addition. All these operations are done approximately with some accuracy $\varepsilon$ . To solve the linear system with the Schur complement and inexact matrix-by-vector products, we used the inexact GMRES as implemented in [5], which implements the adaptive accuracy for the intermediate matrix-by-vector products. Despite the fact that our matrix is positive semidefinite, we decided to use GMRES over the conjugate gradients in this case, since the GMRES method proved to be much more robust to the approximate inversion of the Laplacian.

6 Numerical experiments

As a numerical test, we have selected the following analytical solution:

\begin{split}&u=\frac{1}{4\pi^{2}}\sin 2\pi x(1-\cos 2\pi y),\\ &v=\frac{1}{4\pi^{2}}(1-\cos 2\pi x)\sin 2\pi y,\\ &p=\frac{1}{\pi}\sin 2\pi x\sin 2\pi y.\end{split}

(12)

The corresponding right-hand side is then

\begin{split}&f_{u}=\sin 2\pi y,\\ &f_{v}=\sin 2\pi x,\\ &g=-\frac{1}{\pi}\sin 2\pi x\sin 2\pi y.\end{split}

(13)

The solution for $p$ has a perfect structure (rank-1). However, during the iterative process the intermediate ranks may grow quite a lot, and that increases the complexity. The Table 1 provides the dependence of the total computation time for full low-rank method and for the same method in the full format. The third column gives an error of the approximation with respect to the true solution. For all computations the same threshold $\varepsilon=5\cdot 10^{-9}$ was used that is safe to match the discretization error. It is clear, that the low-rank method is slower for small mode sizes (until $n\leq 256$ ) and is becoming faster due to better scaling in $n$ . The almost linear scaling in $n$ is also observed.

Table 1: Numerical results for the two-dimensional sine example

n	Time (LR/full)	Rel. error in $p$ (LR/full)
64	0.17/0.024	1.2e-03/1.2e-03
128	0.29/0.061	2.9e-04/2.9e-04
256	0.54/0.39	7.8e-05/7.8e-05
512	1.2/2.3	1.9e-05/1.9e-05
1024	2.8/10	8.6e-06/8.2e-06

In Figure 2 we present the approximate rank of the Krylov vectors in the inexact GMRES and the residual computed in the GMRES method. We see that the rank decreases with the iteration number (which is consistent with the inexact GMRES theory, the next Krylov vectors are approximated with less accuracy), and the GMRES converges linearly.

Figure 2: The ranks of Krylov vectors and the convergence of the Uzawa-GMRES method

6.1 Lid-driven cavity flow

The second example is the classical lid-driven cavity test (see the picture). The homogenius Dirichlet boundary conditions are introduced on all the boundaries for $u$ and $v$ except the upper boundary which is moving with constant velocity (see Figure 3).

Figure 3: The lid-driven cavity flow

The results are presented in Table 2. Since now we do not the analytic solution, we can compare the solution obtained from the full format with the solution obtained the the low-rank method. The sublinear scaling in the number of unknowns is clearly visible.

Table 2: Numerical results for the two-dimensional lid-driven cavity flow

n	Time (LR/full)	Rel. error
256	0.89/0.63	6.7e-07
512	2.5/3.2	5.6e-07
1024	6.3/13	5.9e-07

Finally, we will compare the new method with the exponential sums approach. We have taken $n=1024$ and $\varepsilon=1e-7$ . That gave the rank of the inverse $35$ . The typical rank of the Krylov vector for the lid-driven cavity test was around $30$ as well. The timing for the Hadamard product and then rounding was 1.6 seconds, and for the cross method the time was 0.07 seconds, i.e, more than 20 times faster.

7 Conclusions

We have presented a new way to solve Poisson equation in the low-rank tensor format using cross approximation in the frequency domain. The idea is quite simple, but is very effective. Using the new solver we have implemented the first low-rank solver for the two-dimensional Stokes problem. The low-rank Stokes solver is based on the Uzawa-GMRES approach, and is effective and robust in the experiments. In the numerical experiments only a two-dimensional case was considered, since it is the “worst case” for the low-rank methods due to typically high complexity with respect to the ranks. Even in this setting, the low-rank solver outperforms the full-format finite-difference solver starting from $n=512$ . We plan to implement a three-dimensional variant, and the complexity reduction in three dimensions is expected to be much more visible.

References

[1] M. Bebendorf, Approximation of boundary element matrices, Numer. Mathem., 86 (2000), pp. 565–589.
[2] M. Benzi, G. H. Golub, and J. Liesen, Numerical solution of saddle point problems, Acta numerica, 14 (2005), pp. 1–137.
[3] G. Beylkin and M. J. Mohlenkamp, Numerical operator calculus in higher dimensions, Proc. Nat. Acad. Sci. USA, 99 (2002), pp. 10246–10251.
[4] L. de Lathauwer, B. de Moor, and J. Vandewalle, A multilinear singular value decomposition, SIAM J. Matrix Anal. Appl., 21 (2000), pp. 1253–1278.
[5] S. V. Dolgov, TT-GMRES: solution to a linear system in the structured tensor format, Russ. J. Numer. Anal. Math. Modelling, 28 (2013), pp. 149–172.
[6] I. P. Gavrilyuk, W. Hackbusch, and B. N. Khoromskij, Tensor-product approximation to the inverse and related operators in high-dimensional elliptic problems, Computing, (2005), pp. 131–157.
[7] S. A. Goreinov, I. V. Oseledets, D. V. Savostyanov, E. E. Tyrtyshnikov, and N. L. Zamarashkin, How to find a good submatrix, in Matrix Methods: Theory, Algorithms, Applications, V. Olshevsky and E. Tyrtyshnikov, eds., World Scientific, Hackensack, NY, 2010, pp. 247–256.
[8] S. A. Goreinov and E. E. Tyrtyshnikov, The maximal-volume concept in approximation by low-rank matrices, Contemporary Mathematics, 208 (2001), pp. 47–51.
[9] L. Grasedyck, Existence and computation of low Kronecker-rank approximations for large systems in tensor product structure, Computing, 72 (2004), pp. 247–265.
[10] W. Hackbusch and D. Braess, Approximation of $\frac{1}{x}$ by exponential sums in $[1,\infty]$ , IMA J. Numer. Anal., 25 (2005), pp. 685–697.
[11] W. Hackbusch, B. N. Khoromskij, and E. E. Tyrtyshnikov, Approximate iterations for structured matrices, Numer. Mathematik, 109 (2008), pp. 365–383.
[12] B. N. Khoromskij, Structured rank- $(r_{1},\ldots,r_{d})$ decomposition of function-related operators in $\mathbb{R}^{d}$ , Comput. Meth. Appl. Math, 6 (2006), pp. 194–220.
[13] B. N. Khoromskij and V. Khoromskaia, Multigrid accelerated tensor approximation of function related multidimensional arrays, SIAM J. Sci. Comput., 31 (2009), pp. 3002–3026.
[14] E. A. Muravleva, On the kernel of discrete gradient operator, Numer. Meth. Prog., 9 (2008), pp. 93–100.
[15] , Numerical methods based on variational inequalities for viscoplastic Bingham media, PhD thesis, INM RAS, Moscow, 2010.
[16] I. V. Oseledets and E. A. Muravleva, Fast orthogonalization to the kernel of discrete gradient operator with application to the Stokes problem, Linear Algebra Appl., 432 (2010), pp. 1492–1500.
[17] I. V. Oseledets, D. V. Savostianov, and E. E. Tyrtyshnikov, Tucker dimensionality reduction of three-dimensional arrays in linear time, SIAM J. Matrix Anal. Appl., 30 (2008), pp. 939–956.
[18] I. V. Oseledets, D. V. Savostyanov, and E. E. Tyrtyshnikov, Linear algebra for tensor problems, Computing, 85 (2009), pp. 169–188.
[19] I. V. Oseledets and E. E. Tyrtyshnikov, Approximate inversion of matrices in the process of solving a hypersingular integral equation, Comp. Math. and Math. Phys., 45 (2005), pp. 302–313.
[20] , TT-cross approximation for multidimensional arrays, Linear Algebra Appl., 432 (2010), pp. 70–88.
[21] V. Rokhlin and N. Yarvin, Generalized Gaussian quadratures in singular value decompositions of integral operators, SIAM J. Sci. Comput., 20 (1998), pp. 699–718.
[22] L. R. Tucker, Some mathematical notes on three-mode factor analysis, Psychometrika, 31 (1966), pp. 279–311.
[23] E. E. Tyrtyshnikov, Incomplete cross approximation in the mosaic–skeleton method, Computing, 64 (2000), pp. 367–380.