Have diagonalization, need orthogonal diagonalization

Question

I have a symmetric matrix $A$ and I want to find an orthogonal diagonalization of $A$, i.e, matrices $B$ and $D$ such that:

$$ A = B D B^{-1}$$

and $D$ is diagonal and $B^{-1} = B^T$.

I use a math software _(sympy) that lets my diagonalize any matrix, i.e, it lets me find $C$ and $D$ such that:

$$ A = C D C^{-1}$$

and $D$ is diagonal, but usually $C^{-1}\neq C^T$.

Can I use this function to find an orthogonal diagonalization of a symmetric matrix?

Answer 1

I don't know sympy and so I don't know whether or in which form it has a diagonalization by an explicite SVD-decomposition. So here a somehow "pseudocode" how to arrive at that.

---

If $B$ is orthogonal, then $B^{-1}=B^\tau$ . Being orthogonal means $B$ is a rotation-matrix. So if you do a rotation on the rows and the same rotation, but transposed, on the columns then you arrive at a diagonal matrix $D$ and a suitable matrix $B$ .

I have implemented such a procedure as standard-diagonalization procedure for symmetric matrices (main target: correlation-matrices in my statistics-projects) in the matrix-language-program "MatMate". Here, when we have a symmetric matrix we can simply write

A = randomn(4,4)
A = A * A'  // A is now a symmetric matrix
// ---------------------------
B = gettrans(A,"pca")    // rotating a symmetrix matrix to PC-position
                         // getting back the transformation-/rotation matrix
B_I = inv(B)
D = B_I * A * B 
// check:
CHK = B * D * B_1  - A   // to be zero if diagonalization ok
CHK = B_I - B'            // to be zero if B = B^T

$$ \text{ A } = \small \left[\begin{array} {rrrr} 7.2591& -3.5564& -0.5114& 0.0493\\ -3.5564& 3.6224& -2.2029& -0.1632\\ -0.5114& -2.2029& 3.3816& -0.3112\\ 0.0493& -0.1632& -0.3112& 1.7369 \end{array} \right] $$

---

$$ \text{ B }=\small \left[\begin{array} {rrrr} 0.8280& 0.3880& 0.0995& 0.3923\\ -0.5461& 0.3984& 0.1362& 0.7242\\ 0.1263& -0.8280& 0.0273& 0.5456\\ 0.0116& 0.0713& -0.9853& 0.1548 \end{array}\right] $$

$$ \text{ D } \small = \left[ \begin{array} {rrrr} 9.5274& 0.0000& 0.0000& 0.0000\\ 0.0000& 4.7079& 0.0000& 0.0000\\ 0.0000& 0.0000& 1.7631& 0.0000\\ 0.0000& 0.0000& 0.0000& 0.0016 \end{array} \right] $$

$$ \text{ B}^{-1} = \small \left[ \begin{array} {rrrr} 0.8280& -0.5461& 0.1263& 0.0116\\ 0.3880& 0.3984& -0.8280& 0.0713\\ 0.0995& 0.1362& 0.0273& -0.9853\\ 0.3923& 0.7242& 0.5456& 0.1548 \end{array} \right] $$

If you can't reingeneer the gettrans(B,"pca")-procedure (it is also a SVD-decomposition) then a simpler (?) method is the rotation to triangular shape, but this must then be iterated:

D=A
B=einh(4)   \\ initialize the iteratively converging rotation-process
\\ ----------------------
t = gettrans(D,"drei")    // iterate this three steps until convergence
B=B * t                   //
D = t' * D' * t           //  
\\ ----------------------
B_I = inv(B)

This routine gave the same results to 10 decimal digits precision by ca 20 iterations.