LinearSolve on a singular matrix

Question

I have some singular transition rate matrices $m$ (columns add to zero). They have exactly one zero eigenvalue that I want to find the corresponding eigenvector of (the rest of the eigenvalues are negative). This corresponds to the stationary distribution of a continuous-time Markov process. For my actual problem, which is huge and sparse, the fastest method I've found uses (abuses?) LinearSolve. What I don't understand, is why this method works at all.

Here's a minimal example:

m = {{-0.1, 2., 0}, {0.1, -4.1, 6.}, {0, 2.1, -6.}};
Eigenvalues[m]
(* {-8.73318, -1.46682, 0.} *)

Here are two (slower) ways to get the answer I'm looking for:

Eigenvectors[m][[-1]]
(* {0.9986, 0.04993, 0.0174755} *)
NullSpace[m]
(* {{0.9986, 0.04993, 0.0174755}} *)

At first, I thought I could use LinearSolve[m, {0, 0, 0}] but that only gives the trivial solution {0., 0., 0.}. On a lark, I tried replacing the {0, 0, 0} with a small vector. Surprisingly, this works (although with the warning LinearSolve::luc):

x = LinearSolve[m, {10^-10, 10^-10, 10^-10}];
x/Norm[x]
(* {0.9986, 0.04993, 0.0174755} *)

Even more surprising to me, replacing the small vector with a decidedly non-small vector still works with a warning:

x = LinearSolve[m, {1, 1, 1}];
x/Norm[x]
(* {0.9986, 0.04993, 0.0174755} *)

For the record, this answer doesn't actually solve the problem given to LinearSolve (although it solves MY problem):

m.x - {1, 1, 1}
(* {1., -1., -1.} *)

Finally, my original problem with a SparseArray doesn't even generate a warning. This can be reproduced as:

x = LinearSolve[SparseArray[m], {1, 1, 1}];
x/Norm[x]
(* {0.9986, 0.04993, 0.0174755} *)

So, my questions are: why does this work? and, can I count on this working in general?

Answer 1

The example with one zero eigenvector can be explained like this. Of a matrix $M$, let $e_k$ be the unit eigenvectors with nonzero eigenvalues $\lambda_k$ and $n$ be the unit eigenvector with eigenvalue zero. Suppose when we apply LinearSolve to the problem $Mx=b$, we get an approximate solution $x^*$. Suppose $x^*$ is the exact solution to a closely related problem $M^* x = b$, where $M^*$ is a matrix approximately equal to $M$. Let $e^*_k \approx e_k$, $n^* \approx n$ be the unit eigenvectors of $M^*$ with eigenvalues $\lambda^*_k \approx \lambda_k$, $\nu^* \approx 0$. If $$b = c_0 n^* + \sum_{k} c_k e^*_k\,,$$ then the solution $x^*=(M^*)^{-1}b$ will be $$x^* = (\nu^*)^{-1} c_0 n^* + \sum_{k} (\lambda^*_k)^{-1} c_k e^*_k\,.$$ The value of $(\nu^*)^{-1}$ for the OP's m is around $10^{15}$ and much, much bigger than the $(\lambda^*_k)^{-1}$, so that $${x^* \over \|x^*\|} \approx n^* \approx n \quad\text{(up to sign)}\,,$$ unless $c_0$ is quite small.

So the procedure will give a good approximation as long as $c_0$ is not too much smaller than the other coefficients $c_k$ relative to the other eigenvalues. I haven't thought of a way to verify that hypothesis, except to say that it is unlikely to pick a $b$ for which $c_0$ is very small (unless the problem has symmetry and you pick a corresponding $b$). The main issue is whether any of the $\lambda_k$ can be very small. Perhaps that can be known from the particular kind of $M$ that arises.

Answer 2

I am not sure whether yours is a relyable method. But the idea behind it is intriguing. So I tried to derive a method for computing the null space of a matrix from it, but without division by zero. Hopefully, this will be more robust.

Background

The backend of LinearSolve is LAPACK and computes a (permuted) LU-decomposition. The default method ("Multifrontal") for sparse matrices is UMFPACK. It also computes essentially a permuted LU-decomposition, but it uses a clever permutation to make the factors as sparse quite sparse. (I don't say "as sparse as possible" because it's not true; one had to solve a NP-hard problem for that.)

With the input $n\times n$-matrix $A$, a permuted LU-decomposition consists of

• a lower triangular square matrix $L$ with only 1 on the diagonal
• an upper triangular square matrix $U$
• and a permutation of $\{1,\dotsc,n\}$

such that (in Mathematica notation)

 L.U == A[[p]]

The $L$-factor is always invertible. Thus, all information about singularity is contained in the $U$-factor. It is upper triangular, so the rank of $A$ (equalling the rank of $U$) equals the number of nonzero diagonal elements of $U$. I guess that an estimator for the condition number of $A$ is computed somehow from the diagonal elements of $U$ (in a comparatively cheap way compared to computing the minimal and maximal absolute eigenvalue). Anyways, zeros (or numbers close to zero) on the diagonal of $U$ are strong indicators for a singular matrix. Moreover, when solcing linear systems involving U (the first step of solving $L \,U \, x = b$) via back substitution, one has to divide by the diagonal elements. If some of these are close to zero, they shadow somewhat the dependence on the right hand side.

Some alternate algorithm

We can use this knowledge in order to avoid the division by (nearly) zeros. The idea is to modify the factor U to a matrix V that is invertible and to solve linear systems involving V in order to compute a basis of the null space.

QuickNullSpace[A_?SquareMatrixQ, threshold_: 1. 10^-10] := 
  QuickNullSpace[Quiet[LinearSolve[A]], threshold];

QuickNullSpace[S_LinearSolveFunction, threshold_: 1. 10^-10] := 
 Module[{U, L, idx, V, p, q, x},
  U = S["getU"];
  L = S["getL"];
  p = S["getPermutations"];
  {p, q} = If[MissingQ[p], {All, All}, p];
  idx = Flatten[Position[Threshold[Abs[Normal[Diagonal[U]]], threshold], 0.]];
  V = U;
  Do[V[[i, i]] = 1., {i, idx}];
  Transpose[
   Quiet[LinearSolve[
      V,
      SparseArray[
       Transpose[{idx, Range[Length[idx]]}] -> 1., {Length[U], 
        Length[idx]}, 0.]
      ]][[q]]
   ]
  ]

Dense array usage example

Experimentally, this performs quite well for dense matrices. Here we test it on some random $1000 \times 1000$-matrix with $30$-dimensional null space. Some more tweaking is necessary in order to cope with the column permutations that might be needed in order to get sparse factors.

n = 1000;
m = 30;
U = RandomVariate[CircularRealMatrixDistribution[n]];
V = RandomVariate[CircularRealMatrixDistribution[n]];
A = U.(Join[ConstantArray[0., {m}], RandomReal[{-1, 1}, {n - m}]] V);

r1 = QuickNullSpace[A]; // RepeatedTiming // First
r2 = NullSpace[A]; // RepeatedTiming // First
MatrixRank[r1] == m
Max[Abs[A.Transpose[Orthogonalize@r1]]]
Max[Abs[A.Transpose[r2]]]

0.031

0.25

True

1.50924*10^-13

1.28695*10^-16

Sparse matrix usage example

The method seems to work also for sparse matrices. Here an example where it is applied to the graph Laplacian of some random graph with at least $4$ connected components; each component will add another vector to the null space.

G = GraphDisjointUnion @@ Table[RandomGraph[{100, 600} 10], {4}];
A = With[{A = AdjacencyMatrix[G]},
   N[DiagonalMatrix[SparseArray[Total[A, {2}]]] - A]
   ];
r1 = QuickNullSpace[A]; // RepeatedTiming // First
r2 = NullSpace[A]; // RepeatedTiming // First
MatrixRank[r1] == MatrixRank[r2]
Max[Abs[A.Transpose[Orthogonalize@r1]]]
Max[Abs[A.Transpose[r2]]]

0.27

14.

True

7.61197*10^-15

9.51322*10^-15

That's an over $50$-fold speedup. Moreover, QuickNullSpace will perform relatively better with large matrices of increasing size.

Remarks

• This implementation is highly experimental. No guarantees from my side. Suggestions are welcome. I don't know how reliable (how accurate, how stable) this method is.

• Notice that the QuickNullSpace[A] is not Orthogonalized. I prefer to leave the choice to the user.

  • There is some potential for further speedup. Currently, two calls to LinearSolve are necessary: One for the factorization of A and one for solving several linear systems with the V which we obtains as perturbation of U. Since V is an upper triangular matrix, we could in principle skip the factorization step and start immediately with the back substitution. For the dense matrix case, the solution might lie around somewhere in the context "LinearAlgebra`LAPACK`". But I have not clue how to do it for the sparse matrix case without resorting to LibraryLink calls.