# Solving shifted linear system

**URL:** <https://discourse.julialang.org/t/solving-shifted-linear-system/111265>\
**Category:** Numerics\
**Tags:** linearalgebra\
**Created:** [March 6, 2024, 10:11pm UTC](https://discourse.julialang.org/t/solving-shifted-linear-system/111265 "2024-03-06T22:11:48Z")\
**Posts on this page:** 1\
**Showing post:** 9

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 10, 2024, 6:20pm UTC](https://discourse.julialang.org/t/solving-shifted-linear-system/111265/9 "2024-03-10T18:20:04Z")

</div>

> [@mstewart](#):
>
> Cholesky factorization of T +k^2 I s O(n)

You don’t even need an explicit Cholesky factorization, since for symmetric tridiagonal matrices I seem to recall that we have a `T \ b` method that works directly with `T` (and is also O(n)). In fact, there is an `ldiv!(T, b; shift=µ)` method that computes `(T + µ*I) \ b` in-place in `b`. _Update:_ [checking the source code](https://github.com/JuliaLang/julia/blob/5c7d24493ebab184d4517ca556314524f4fcb47f/stdlib/LinearAlgebra/src/tridiag.jl#L250) it seems that it does form an `ldlt` factorization explicitly, however, though it is still true that it is O(n); in principle this could probably be elided into the solve.

If you compute a Hessenberg factorization object `F = hessenberg(A)` in Julia, this is exploited under the hood. As [explained in the documentation](https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.Hessenberg) it directly supports shifts: `(F + µ*I) \ b` or `ldiv!(F + µ*I, b)` is very cheap (O(n) for Hermitian matrices, and O(n^2) for general matrices).

What is (currently) missing is a BandedMatrices.jl method that creates a `Hessenberg` factorization object from a banded matrix, with an implicit representation of Q as a product of Givens rotations. This should be straightforward to add, in principle, since much of the code is already written.

The Hessenberg approach could be an especially big win if you are solving for different shifts with the _same right-hand side_, since in that case you may only need to apply the Q transformation _once_ and then just do tridiagonal solves for each shift.

> [@rvaradhan](#):
>
> I know the spectral and SVD decompositions where there is the similar idea that the diagonal shift is absorbed in the middle diagonal matrix

Finding the Hessenberg factorization is significantly cheaper than computing the diagonalization of the matrix — in fact, the first step of diagonalization is usually a Hessenberg factorization, so this saves you the subsequent iterative diagonalization steps.

(For general square matrices, the analogous thing to diagonalization for shifted solves would typically be a Schur factorization, not an SVD, but again Hessenberg is the first step of a Schur calculation so working with Hessenberg is cheaper.)

But again, for banded matrices with small bandwidth, _every_ reasonable algorithm is linear O(n) time, including just re-computing the Cholesky factorization for each shift — so it boils down to a battle of the constant factors and which method is more optimized. It’s not obvious _a priori_ which method will “win”, but on the other hand they are all so fast that it may not matter much.

---

_[View the full topic](https://discourse.julialang.org/t/solving-shifted-linear-system/111265)._