Suppose A
is a complex n x n
matrix where n
is quite big, given as a numeric array. For many values of the complex number z
, and many complex vectors b
with n
components, I need to compute LinearSolve[A-z*IdentityMatrix[n],b]
. But calling LinearSolve
repeatedly is slow, and I need something faster. Note that the matrix Inverse[A-z*IdentityMatrix[n]]
is also known as the resolvent of A
.
In my application, A
is not diagonalizable, so I cannot make an eigendecomposition. But I can always assume that z
is not an eigenvalue of A
. I think a reasonable algorithm is to use a Schur decomposition as a pre-processing step, and to perform a backsubstitution each time the resolvent is applied:
(* random example *)
n = 2000;
A = RandomComplex[{-1-I,1+I},{n,n}];
(* slow implementation *)
resolventSlow[z_,b_] := LinearSolve[A-z*IdentityMatrix[n],b];
(* preprocessing *)
{Q,T} = SchurDecomposition[A,RealBlockDiagonalForm->False];
(* faster implementation *)
resolventFast[z_,b_] := Q.LUBackSubstitution[{T-z*IdentityMatrix[n],Range[1,n],1.},Conjugate[Conjugate[b].Q]];
(* timing *)
z = RandomComplex[{-1-I,1+I}];
b = RandomComplex[{-1-I,1+I},n];
RepeatedTiming[x1 = resolventSlow[z,b];] (* about 0.5 seconds *)
RepeatedTiming[x2 = resolventFast[z,b];] (* about 0.1 seconds *)
Chop[Norm[x1-x2]] (* zero as expected *)
The function LUBackSubstitution
that I use for backsubstitution has been obsolete since 2003, and its specification is a little arcane. See for example this v4 documentation and this documentation of LUDecomposition.
Question: Is there an idiomatic replacement for LUBackSubstitution
in the code above, that does not explicitly invoke obsolete symbols or low-level libraries such as BLAS? Alternatively, is there another useful algorithm for the resolvent that I could use instead?
LinearSolveFunction
usingLinearSolve[T-z*IdentityMatrix[n]]
sinceT
is already upper triangular. $\endgroup$LUBackSubstitution
is going anywhere so it should be safe to use, deprecation notwithstanding. $\endgroup${qq, tt} = SchurDecomposition[A, RealBlockDiagonalForm -> False]; resolventFast2[z_, b_] := Module[{tr = tt - z IdentityMatrix[n], bb = b . Conjugate[qq]}, LinearAlgebra`BLAS`TRSV["U", "N", "N", tr, bb]; qq . bb]
$\endgroup$TRSV
code works for me and is about as fast as theLUBackSubstitution
code. I explicitly excluded BLAS since I was hoping there could be, for example, some way to giveLinearSolve
a hint about upper triangularity. $\endgroup$UpperTriangularMatrix[]
, which you can try for your resolvent computation. $\endgroup$