Symmetric positive-definite preconditioning (source code)

= Symmetric positive-definite preconditioning
{title2=$B=M^{-1/2}AM^{-1/2},\quad y=M^{1/2}x$}

Let $A$ and $M$ be real <symmetric positive-definite matrices>. The <principal square root of a positive semidefinite matrix> construction gives an <invertible matrix> for each. Set
$$
S=M^{1/2},\qquad B=S^{-1}AS^{-1},\qquad y=Sx,\qquad c=S^{-1}b.
$$
Then $Ax=b$ is equivalent to $By=c$. This is a <matrix congruence>, and $B$ is a <symmetric positive-definite matrix>: $B^T=B$, and for every nonzero <vector> $v$,
$$
v^TBv=(S^{-1}v)^TA(S^{-1}v)>0,
$$
since $S^{-1}v\ne0$. This short proof establishes both requirements for applying the <Conjugate gradient method> to the transformed <linear system>.

The left transformed <matrix> has the same <eigenvalues> because
$$
M^{-1}A=S^{-1}BS.
$$
Its correct geometry uses the <weighted inner product> $\langle u,v\rangle_M=u^TMv$: $\langle u,M^{-1}Av\rangle_M=u^TAv=\langle M^{-1}Au,v\rangle_M$. Thus it is symmetric in this <inner product>, even when it is not a <symmetric matrix> in ordinary coordinates. The spectral <condition number> governing the transformed symmetric problem is $\kappa_2(B)$, not generally $\kappa_2(M^{-1}A)$.

The <Preconditioned conjugate gradient method> realizes this transformation using auxiliary solves $Mz_k=r_k$, without constructing $S$ or its <matrix inverse>. With fixed $M$, initialize $r_0=b-Ax_0$ and $p_0=z_0$; while $r_k\ne0$, use
$$
\alpha_k=\frac{r_k^Tz_k}{p_k^TAp_k},\quad x_{k+1}=x_k+\alpha_kp_k,\quad r_{k+1}=r_k-\alpha_kAp_k,
$$
then solve $Mz_{k+1}=r_{k+1}$ and set
$$
\beta_k=\frac{r_{k+1}^Tz_{k+1}}{r_k^Tz_k},\qquad p_{k+1}=z_{k+1}+\beta_kp_k.
$$
Stop at a zero <linear-system residual> before forming another ratio. Positive definiteness makes the denominators positive before termination.<Math 221, Matrix Computations, Fall 2023/Notes for Ma221 Lecture 15, Dec 1, 2023>[Math 221, Matrix Computations, Fall 2023; Notes for Ma221 Lecture 15, Dec 1, 2023; PDF p. 1, equivalent symmetric system and Preconditioned CG algorithm]{ref}

An <invertible matrix> alone is insufficient: with $A=I$ and $M=\operatorname{diag}(-1,1)$, the left transformed <matrix> has <eigenvalues> $-1,1$, and $M$ has no real <symmetric positive-definite matrix> square root. Even a valid <symmetric positive-definite matrix> can worsen conditioning if its scaling is unsuitable, as the next worked example shows.