Symmetric LU

We begin with a symmetric \(A\).

A = [  2     4     4     2
       4     5     8    -5
       4     8     6     2
       2    -5     2   -26 ];

Carrying out our usual elimination in the first column leads us to

using LinearAlgebra
L1 = diagm(0=>ones(4))
L1[2:4,1] = [-2,-2,-1]
A1 = L1*A
4×4 Array{Float64,2}:
 2.0   4.0   4.0    2.0
 0.0  -3.0   0.0   -9.0
 0.0   0.0  -2.0   -2.0
 0.0  -9.0  -2.0  -28.0

But now let’s note that if we transpose this result, we have the same first column as before! So we could apply again and then transpose back.

A2 = (L1*A1')'
4×4 Adjoint{Float64,Array{Float64,2}}:
 2.0   0.0   0.0    0.0
 0.0  -3.0   0.0   -9.0
 0.0   0.0  -2.0   -2.0
 0.0  -9.0  -2.0  -28.0

Using transpose identities, this is just

A2 = A1*L1'
4×4 Array{Float64,2}:
 2.0   0.0   0.0    0.0
 0.0  -3.0   0.0   -9.0
 0.0   0.0  -2.0   -2.0
 0.0  -9.0  -2.0  -28.0

Now you can see how we proceed down and to the right, eliminating in a column and then symmetrically in the corresponding row.

L2 = diagm(0=>ones(4))
L2[3:4,2] = [0,-3]
A3 = L2*A2*L2'
4×4 Array{Float64,2}:
 2.0   0.0   0.0   0.0
 0.0  -3.0   0.0   0.0
 0.0   0.0  -2.0  -2.0
 0.0   0.0  -2.0  -1.0

Finally, we arrive at a diagonal matrix.

L3 = diagm(0=>ones(4))
L3[4,3] = -1
D = L3*A3*L3'
4×4 Array{Float64,2}:
 2.0   0.0   0.0  0.0
 0.0  -3.0   0.0  0.0
 0.0   0.0  -2.0  0.0
 0.0   0.0   0.0  1.0