Algorithm for the resolvent of a matrix, and LUBackSubstitution

Question

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?

Answer 1

As pointed out by @J.M., starting with Version 13.1 (June 2022) there is UpperTriangularMatrix. It is efficient with LinearSolve, see the Details-section of the documentation.

Therefore one can now use the following idiomatic code:

(* V13.1 *)

(* preprocessing *)
{Q,T} = SchurDecomposition[A,RealBlockDiagonalForm->False];

(* new solution *)
resolventFast[z_,b_] := Q.LinearSolve[UpperTriangularMatrix[T-z*IdentityMatrix[n]],
                                      Conjugate[Conjugate[b].Q]];

Possible improvement 1. To also exploit efficient storage, one could call T = UpperTriangularMatrix[T] right after SchurDecomposition.

Possible improvement 2. Version 13.1 also has BlockUpperTriangularMatrix (experimental). If A is real then one can use SchurDecomposition in default mode which produces a block-upper-triangular T, but has the advantage that both Q and T are real which saves memory:

(* V13.1 -- only for real A *)

{Q,T} = SchurDecomposition[A];
resolventFast[z_,b_] := Q.LinearSolve[BlockUpperTriangularMatrix[
                                        T-z*IdentityMatrix[n]],b.Q];