4
$\begingroup$

Let $A$ be a general symmetric matrix operator and $P$ be the unique orthogonal projection onto $\operatorname{Range}(A) = \operatorname{Null}(A)^\perp$.

Analytically, the system $$Ax = Pb$$ should have a solution for any vector $b$. If we further stipulate that this solution lie in $\operatorname{Range}(A)$, the solution should be unique. (And in the positive definite case, this is precisely the solution that the conjugate gradient algorithm will find in exact arithmetic.)

The above is a common tactic when solving the Poisson equation with either Neumann or periodic boundary conditions; if $A$ is taken to be the Neumann or periodic Laplacian, then $$ P = I - \frac{1}{N}(1_N)(1_N)^T,$$ where $1_N = [1, \cdots, 1]^T$, and we can solve $Ax = Pb$ using a conjugate gradient approach.

But what if the projection operator $P$ doesn't have a simple analytic expression? Can it be calculated numerically from $A$? Is there a way to solve this type of singular problem in the case that $A$ is a black-box symmetric operator?

I realize one approach would be to solve the normal equations $A^2x = Ab$, which avoids using a projection operator entirely, but this also severely worsens the conditioning of the system, so I'm interested to know if there's another way.

$\endgroup$
3
$\begingroup$

Computing the null space without any further knowledge is expensive in general. There are some iterative methods that converge to minimum-norm solutions even when presented with inconsistent right hand sides. Choi, Paige, and Saunders' MINRES-QLP is a nice example of such a method. For non-symmetric problems, see Reichel and Ye's Breakdown-free GMRES.

In practice, usually some characterization of the null space is important for effective preconditioning. Since most practical problems require preconditioning, the purely iterative methods have seen limited adoption. Note that in case of very large null space, preconditioners will often be used in an auxiliary space where the null space has been removed. See the "auxiliary-space Maxwell" methods for more details.

$\endgroup$
2
$\begingroup$

Jed has already answered your question, implicitly assuming that your matrix is large and sparse. If this should not be the case, you can do an eigenvalue decomposition and write $A=X\Lambda X^T$ where the columns of $X$ correspond to the eigenvectors and the entries of the diagonal matrix $\Lambda$ are the eigenvalues. (Note that this decomposition is not easily computable if $A$ is large, as Jed mentioned.) Then, the null space consists of exactly those vectors spanned by the eigenvectors corresponding to zero (or very small) eigenvalues, i.e., the projector onto the null space can be written as $$ 1_N = X \tilde\Lambda X^T $$ where $(\tilde\Lambda)_{ii} = 0$ if $|\Lambda_{ii}|\ge \varepsilon$ and $(\tilde\Lambda)_{ii} = \Lambda_{ii}$ if $|\Lambda_{ii}|<\varepsilon$.

$\endgroup$
  • 1
    $\begingroup$ This procedure will work. In practice, an SVD or QR decomposition is typically preferred, in which case left singular vectors corresponding to sufficiently small singular values make up a basis for the nullspace. $\endgroup$ – Geoff Oxberry Jun 2 '13 at 6:24
  • 1
    $\begingroup$ @GeoffOxberry: Yes, you are correct. I meant to add that (the above is the conceptual view, not how it's implemented) but forgot. That said, if $A$ is a symmetric matrix, isn't the SVD the same as the eigenvalue decomposition, and isn't QR just one way of computing it? $\endgroup$ – Wolfgang Bangerth Jun 2 '13 at 17:30
  • $\begingroup$ Yes. All of the above also apply to the general case, although of course the SVD and eigendecomposition are not necessarily equal if $A$ is not Hermitian. Although SVD can be calculated using QR (and vice versa), I've been told the SVD is the more accurate algorithm. $\endgroup$ – Geoff Oxberry Jun 2 '13 at 22:23
  • 2
    $\begingroup$ To follow up, I looked this up, and the SVD is only the same as the eigenvalue decomposition if $A$ is Hermitian and positive semidefinite. If $A$ is Hermitian, but not positive semidefinite, it still has orthogonal eigenvectors, which can be used to calculate a projector. $\endgroup$ – Geoff Oxberry Jun 4 '13 at 5:15
  • $\begingroup$ Thanks for double checking. I will admit that the many classes of matrices and what that implies for their eigenvalues/eigenvectors is not my forte :-( $\endgroup$ – Wolfgang Bangerth Jun 5 '13 at 2:42

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.