# Solving a Shifted Linear System of 2 Symmetric Positive Semi Definite (SPSD) Matrices

**URL:** <https://discourse.julialang.org/t/solving-a-shifted-linear-system-of-2-symmetric-positive-semi-definite-spsd-matrices/131636>\
**Category:** Optimization (Mathematical)\
**Tags:** linearalgebra, convex-optimization, linear-regression, linearsolve, least-squares\
**Created:** [August 16, 2025, 10:47am UTC](https://discourse.julialang.org/t/solving-a-shifted-linear-system-of-2-symmetric-positive-semi-definite-spsd-matrices/131636 "2025-08-16T10:47:40Z")\
**Posts on this page:** 1\
**Showing post:** 3

<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:** [August 16, 2025, 4:27pm UTC](https://discourse.julialang.org/t/solving-a-shifted-linear-system-of-2-symmetric-positive-semi-definite-spsd-matrices/131636/3 "2025-08-16T16:27:35Z")

</div>

> [@mstewart](#):
>
> What I have done before is to just to give up on exploiting symmetry.

One option is to first Cholesky-factor C = U^T U, let y = U x, then rewrite the system as:

(\underbrace{U^{-T} A U^{-1}}\_{B} + \lambda I) y = U^{-T} A^T b = c

Note that the matrix B is still SPSD. Note also that, as usual, you don’t actually invert any matrices: you compute B via triangular solves `B = U' \ A / U`.

Once you have it in this form, you can call `hessenberg(Hermitian(B))` to factorize B into Q H Q^T, where H is symmetric tridiagonal. The advantage of this is that you can then do shifted solves without repeating the factorization, since B + \lambda I = Q (H + \lambda I) Q^T.

(Diagonalizing B into eigenvectors would also accomplish this, but is more expensive. On the other hand, diagonalization results in a faster solve step, so it depends on how many \lambda you need to solve with.)

Moreover, the Hessenberg factorization objects in Julia have [efficient specialized support for shifted solves](https://github.com/JuliaLang/julia/pull/31853), for exactly this sort of application.

In short, you can do something like:

```julia-auto
# pre-processing
U = cholesky(C).U # UpperTriangular
B = Hermitian(U' \ A / U)
c = U' \ (A'b)
F = hessenberg(B) # tridiagonal

# solve for each λ
λ = 0.1234
y = (F + λ*I) \ c
x = U \ y

# check:
@show x ≈ (A + λ*C) \ (A'b) # should print "true"

```

Note that F + \lambda\*I is actually lazy — it doesn’t even allocate a new matrix. And you can do the y and x solves in-place too, if you have pre-allocated vectors `x` and `y`:

```julia-auto
ldiv!(y, F + λ*I, c)
ldiv!(x, U, y)

```

> [@RoyiAvital](#):
>
> If there a specialized way to solve this for various \lambda values without the squared condition number, it would be even better.

Usually if you are doing Tikhonov regularization, the resulting system should be well-conditioned, no?

If it’s ill-conditioned, you could equivalently solve the ordinary least-squares problem

\min\_x \left\Vert \begin{pmatrix} D \\ \sqrt{\lambda} E \end{pmatrix} x - \begin{pmatrix} b \\ 0 \end{pmatrix} \right\Vert\_2

by QR as usual (`\` in Julia), but I’m not sure if there is a way to re-use calculations for multiple \lambda in that case.

---

_[View the full topic](https://discourse.julialang.org/t/solving-a-shifted-linear-system-of-2-symmetric-positive-semi-definite-spsd-matrices/131636)._