2
$\begingroup$

I have a linear system of equations $Ax = b$ where the number of unknowns $N$ is intractably large but the right-hand side has only small support and the unknown $x$ is known to decay exponentially. I hence want to find a sequence of reduced linear systems $A_n x_n = b_n$ with only $n \leq N$ unknowns such that $x_n \to x$ as $n \to N$. Is there any theory available on how to find suitable $A_n, b_n$ such that I can guarantee well-posedness of the reduced problems and convergence, possibly under additional assumptions on $A$ and $b$?

I have a fairly limited understanding of what wavelets are and how they work, but as far as I see we should face the same problem there: we have an in principle infinite expansion of the solution in terms of the wavelet basis functions, but the smoothness of the solution implies that the corresponding coefficients decay and hence we can truncate. This makes me wonder how the above problems of well-posedness and convergence for truncated systems are addressed there.

$\endgroup$
2
$\begingroup$

Your question is phrased in a very general way and so it is hard to recommend a specific way to go. I will give a recommendation that might help you.

It is well known that every matrix $A\in\mathbf{R}^{n\times n}$ can be factored using the Singular Value Decomposition so that $A=U\Sigma V^T$ where $U$ and $V$ are orthogonal matrices and $\Sigma$ is diagonal. There are standard computational methods for the SVD. Since $U$ and $V$ are orthogonal, the equation $Ax=b$ can be written $x = V\Sigma^{-1}U^T b$. So far, nothing is easier since all we have done is find another way of finding the inverse of $A$. The goal is to approximate $V\Sigma^{-1}U^T$ in an efficient way.

Luckily there are methods to compute the truncated SVD where only the first $k$ columns of $U$ and rows of $V$ are computed. It can be shown that, at least in some sense, that is the best approximation that you can find.

You might search for 'model reduction' or 'model reduction svd' to get more of a view of how the SVD has been used by others to construct a reduced model.

$\endgroup$
1
$\begingroup$

Here are two more precise formulations of the problem which will allow for an approximation scheme as asked for by the OP.


Let $A \in \mathbb{K}^{N \times N}$ be a sparse matrix and well-conditioned matrix and $b \in \mathbb{K}^N$ a sparse vector. Then we can associate $A$ with a graph $G = (V,E)$ where $$ V := \{1, \ldots, N\}, \qquad E := \{(i,j) \in V^2 \mid A(i,j) \neq 0\}, $$ which in turns allows us to define a distance $d(i,j)$ between indices given by the minimal number of edges in $G$ we need to take in order to get from $j$ to $i$ (i.e., $d(i,j)$ is the graph distance in $G$). Then one can show $$ |A^{-1}(i,j)| \leq C \, \exp\big(-\gamma \, d(i,j)\big) \qquad \forall i,j \in V $$ for some $\gamma > 0$, see e.g. "Decay Rates for Inverses of Band Matrices" by Demko, Moss and Smith (1984).

Given the above, it is clear that $x = A^{-1}b$ decays exponentially if we arrange its coefficients according to the graph distance from the non-zero entries of $b$, which is the ordering we assume for the following construction of an approximation $x_n = A_n^{-1} b$ with only $n \ll N$ degrees of freedom. Let $A_n$ be the matrix in which all interactions between entries $\leq n$ and $> n$ are dropped, i.e. $$ A_n(i,j) := \begin{cases} A(i,j) & \text{if } i,j \leq n \text{ or } i,j > n, \\ 0 & \text{otherwise}. \end{cases} $$ It is clear that then $x_n(i) = 0$ for $i > n$, and hence it only remains to show that $x_n(i) \approx x(i)$ for $i \leq n$. This can be seen by expressing the error as $$ x - x_n = (A^{-1} - A_n^{-1})\, b = A^{-1} \, (A_n - A) \, A_n^{-1} \, b. $$ We now assume that $A_n^{-1}$ satisfies the same decay estimate as $A^{-1}$, which boils down to assuming that the conditioning of $A_n$ is not worse than that of $A$. This is trivial if $A$ is positive definite but may be hard to verify otherwise, yet it is an assumption you would probably want to make anyways if you want to numerically solve $A_n x_n = b$. Under this assumption, we estimate the error as $$ |x(i) - x_n(i)| \leq C \, \sum_{k > n} \sum_{\ell \leq n} \sum_{j \in \mathrm{supp}(b)} A(k,\ell) \, \exp\big( - \gamma (d(i,k) + d(\ell,j)) \big), $$ which is at most $C \, Nn|\mathrm{supp}(b)| \, \exp(-\gamma \, d(n+1,\mathrm{supp}(b))\big)$ and can be as small as $C \, Nn|\mathrm{supp}(b)| \, \exp(-2\, \gamma \, d(n+1,\mathrm{supp}(b))\big)$ for $i$ close to $1$.


Alternatively, if $A$ is not sparse but $x$ is known to be decaying for other reasons, we can still order the degrees of freedom such that $|x(i)|$ is decreasing for increasing $i$. We then define the projector onto the first $n$ degrees of freedom,

$$ P_n : \mathbb{K}^N \to \mathbb{K}^N, \qquad a \mapsto \begin{cases} a(i) & \text{if } i \leq n,\\ 0 & \text{otherwise}, \end{cases} $$ and set $A_n := P_n A P_n$ and $x_n := A_n^\dagger b$, where $A_n^\dagger$ is the pseudo-inverse of $A_n$ satisfying $A_n^\dagger A_n = A_n A_n^\dagger = P_n$. Note that once again $x_n$ is formally an element of $\mathbb{K}^N$ but has only $n$ non-zeros. We can then write the error as $x_n - x = (x_n - P_n x) + (P_n x - x)$. The term $P_n x - x$ is easily estimated using the decay of $x$, and for the other term we get \begin{aligned} x_n - P_n x &= P_n \, (A_n^\dagger - A^{-1}) \, b \\&= A_n^\dagger \, (A - A_n) \, A^{-1} \, b \\&= A_n^\dagger \, (A - A_n) \, x \\&= A_n^\dagger \, A \, (x - P_n x) \end{aligned} and can hence be estimated as well by the projection error $x - P_n x$ if $A_n^\dagger$ and $A$ are bounded in suitable norms.

$\endgroup$
0
$\begingroup$

This is akin to having a partial differential equation that is posed to an infinite domain. But we can't solve on infinite domains, so we need to truncate it somewhere, and we can typically do that because the source is localized in a small area of the domain and the solution decays away from the source. One then needs to pose an "artificial boundary condition", for example by only solving in a ball of radius $R$ around the center of the source, and prove that the solution converges as $R\rightarrow \infty$.

Without knowing more about the problem you have, there is little I can offer in terms of concrete suggestions. But you could look into the literature on artificial boundary conditions, truncated domains, and exterior problems, to see how this is done there and how that may apply to your own problem.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.