Symmetric Pair Decomposition
An esoteric fact about matrices is that any real matrix can be written as the product of two symmetric matrices. I've known about this fact for years, but never seriously explored the computational aspects. So I'm using this post to clarify my own understanding of what I'll call the symmetric pair decomposition. It turns out that there are open questions. I don't think we know how to reliably compute the factors. But I also have to admit that, even if we could compute them, I don't know of any practical use.
Contents
Theorem
Not many people know about this theorem.
Theorem: Any real matrix is equal to the product of two real symmetric matrices.
Almost Proof
At first glance this theorem has nothing to do with eigenvalues. But here is the beginning of a proof, and a possible algorithm. Suppose that a real matrix $A$ has real, distinct eigenvalues. Then it can be diagonalized by the matrix $V$ of its eigenvectors.
$$A = V D V^{-1}$$
Because we are assuming there are no multiple eigenvalues, the matrix $V$ exists and is nonsingular, and the matrix $D$ is real and diagonal. Let
$$S_1 = V D V^T$$
$$S_2 = V^{-T} V^{-1}$$
Then $S_1$ and $S_2$ are real, symmetric, and their product is
$$S_1 S_2 = A$$
This argument is not a proof. It just makes the theorem plausible. The challenge comes when the matrix has repeated eigenvalues and lacks a full set of eigenvectors, so it cannot be diagonalized. A complete proof would transform the matrix to its Rational Canonical Form or its Jordan Canonical form and construct explicit symmetric factors for the blocks in the canonical form.
Degrees of Freedom
If $A$ is $n$ -by- $n$ and
$$A = S_1 S_2$$
where each of the symmetric matrices has $n(n+1)/2$ independent elements, then this is $n^2$ nonlinear equations in $n^2+n$ unknowns. It looks like there could be an $n$-parameter family of solutions. In my almost proof, each eigenvector is determined only up to a scale factor. These $n$ scale factors show up in $S_1$ and $S_2$ in complicated, nonlinear ways. I suspect that allowing complex scale factors parameterizes the complete set of solutions, but I'm not sure.
Moler's Rules
My two Golden Rules of computation are:
- The hardest things to compute are things that do not exist.
- The next hardest things to compute are things that are not unique.
For the symmetric pair decomposition, our obscure theorem says the decomposition exists, but the degrees of freedom observation says it is probably not unique. Worse yet, the only algorithm we have requires a full set of eigenvectors, which may not exist. We will have to worry about these things.
LU Decomposition
The most important decomposition in numerical linear algebra, the one we use to solve systems of simultaneous linear equations, is the LU decomposition. It expresses a permuted matrix as the product of two triangular factors.
$$P A = L U$$
The permutation matrix $P$ gives us existence and numerical stability. Putting ones on the diagonal of $L$ eliminates $n$ degrees of freedom and give us uniqueness.
Magic Square
Our first example involves one of my favorite matrices.
A = magic(3)
A =
8 1 6
3 5 7
4 9 2
Use the Symbolic Toolbox to compute the eigenvalues and vectors exactly.
[V,D] = eig(sym(A))
V = [ (2*6^(1/2))/5 - 7/5, - (2*6^(1/2))/5 - 7/5, 1] [ 2/5 - (2*6^(1/2))/5, (2*6^(1/2))/5 + 2/5, 1] [ 1, 1, 1] D = [ -2*6^(1/2), 0, 0] [ 0, 2*6^(1/2), 0] [ 0, 0, 15]
Notice that the elements of V and D involve $\sqrt{6}$ and so are irrational. Now let
S1 = simplify(V*D*V') S2 = simplify(inv(V*V'))
S1 = [ 1047/25, -57/25, 27/5] [ -57/25, 567/25, 123/5] [ 27/5, 123/5, 15] S2 = [ 3/16, 1/12, 1/16] [ 1/12, 1/2, -1/4] [ 1/16, -1/4, 25/48]
The $\sqrt{6}$ has disappeared. You can see that S1 and S2 are symmetric, have rational entries, and, as advertised, their product is
Product = S1*S2
Product = [ 8, 1, 6] [ 3, 5, 7] [ 4, 9, 2]
Let's play with the scale factors a bit. I particularly like
V(:,3) = 2 S1 = simplify(V*D*V')/48 S2 = 48*simplify(inv(V*V'))
V = [ (2*6^(1/2))/5 - 7/5, - (2*6^(1/2))/5 - 7/5, 2] [ 2/5 - (2*6^(1/2))/5, (2*6^(1/2))/5 + 2/5, 2] [ 1, 1, 2] S1 = [ 181/100, 89/100, 21/20] [ 89/100, 141/100, 29/20] [ 21/20, 29/20, 5/4] S2 = [ 5, 0, -1] [ 0, 20, -16] [ -1, -16, 21]
Now S1 has decimal fraction entries, and S2 has integer entries, including two zeros. Let's leave the symbolic world.
S1 = double(S1) S2 = double(S2) Product = S1*S2
S1 =
1.8100 0.8900 1.0500
0.8900 1.4100 1.4500
1.0500 1.4500 1.2500
S2 =
5 0 -1
0 20 -16
-1 -16 21
Product =
8 1 6
3 5 7
4 9 2
I can't promise to get such pretty results with other examples.
Nearly Defective
Suppose I want to compute the symmetric pair decomposition of this perturbation of a Jordan block.
e = sym('e','positive'); A = [2 1 0 0; 0 2 1 0; 0 0 2 1; e 0 0 2]
A = [ 2, 1, 0, 0] [ 0, 2, 1, 0] [ 0, 0, 2, 1] [ e, 0, 0, 2]
Here is the eigenvalue decomposition. The vectors have been scaled so that the last component is equal to 1. The eigenvalues are located on a circle in the complex plane centered at 2, with a radius of e^(1/4), which is the signature of an eigenvalue of multiplicity 4.
[V,D] = eig(A); V = simplify(V) D = simplify(D)
V = [ -i/e^(3/4), i/e^(3/4), 1/e^(3/4), -1/e^(3/4)] [ -1/e^(1/2), -1/e^(1/2), 1/e^(1/2), 1/e^(1/2)] [ i/e^(1/4), -i/e^(1/4), 1/e^(1/4), -1/e^(1/4)] [ 1, 1, 1, 1] D = [ 2 - e^(1/4)*i, 0, 0, 0] [ 0, e^(1/4)*i + 2, 0, 0] [ 0, 0, e^(1/4) + 2, 0] [ 0, 0, 0, 2 - e^(1/4)]
Here is the symmetric pair decomposition resulting from this eigenvalue decomposition.
S1 = simplify(V*D*V.') S2 = simplify(inv(V*V.'))
S1 = [ 0, 4/e, 8/e, 0] [ 4/e, 8/e, 0, 0] [ 8/e, 0, 0, 4] [ 0, 0, 4, 8] S2 = [ 0, 0, e/4, 0] [ 0, e/4, 0, 0] [ e/4, 0, 0, 0] [ 0, 0, 0, 1/4]
Well, this sort of does the job. S1 and S2 are symmetric and their product is equal to A.
Product = S1*S2
Product = [ 2, 1, 0, 0] [ 0, 2, 1, 0] [ 0, 0, 2, 1] [ e, 0, 0, 2]
But I am worried that the factors are very badly scaled. As I make e smaller, the large elements in S1 get larger, and the small elements in S2 get smaller. The decomposition breaks down.
A Better Decomposition
A better decomposition is just a rotation. These two matrices are symmetric and their product is A.
S2 = sym(rot90(eye(size(A)))) S1 = A/S2 Product = S1*S2
S2 = [ 0, 0, 0, 1] [ 0, 0, 1, 0] [ 0, 1, 0, 0] [ 1, 0, 0, 0] S1 = [ 0, 0, 1, 2] [ 0, 1, 2, 0] [ 1, 2, 0, 0] [ 2, 0, 0, e] Product = [ 2, 1, 0, 0] [ 0, 2, 1, 0] [ 0, 0, 2, 1] [ e, 0, 0, 2]
Can I reproduce this decomposition by rescaling the eigenvectors? Here is code that uses the symbolic solve function to compute new scale factors. If you want to see how it works, download this M-file using the link at the end of this post, remove the semicolons in this section, and run or publish it again.
s = sym('s',[4,1]);
V = V*diag(s);
T = simplify(inv(V*V.'));
soln = solve(T(:,1)-S2(:,1));
s = [soln.s1(1); soln.s2(1); soln.s3(1); soln.s4(1)]
s =
(e^(3/4)*i)^(1/2)/2
(-e^(3/4)*i)^(1/2)/2
e^(3/8)/2
(-e^(3/4))^(1/2)/2
These scale factors are complex numbers with magnitude e^(3/8)/2. Let's rescale the eigenvectors. Of course, the eigenvalues don't change.
[V,D] = eig(A); V = simplify(V*diag(s))
V = [ -(-1)^(3/4)/(2*e^(3/8)), (-1)^(1/4)/(2*e^(3/8)), 1/(2*e^(3/8)), -i/(2*e^(3/8))] [ -(-1)^(1/4)/(2*e^(1/8)), -1/(-1)^(1/4)/(2*e^(1/8)), 1/(2*e^(1/8)), i/(2*e^(1/8))] [ ((-1)^(3/4)*e^(1/8))/2, -((-1)^(1/4)*e^(1/8))/2, e^(1/8)/2, -(e^(1/8)*i)/2] [ ((-1)^(1/4)*e^(3/8))/2, (1/(-1)^(1/4)*e^(3/8))/2, e^(3/8)/2, (e^(3/8)*i)/2]
Now these eigenvectors produce the same stable decomposition as the rotation.
S1 = simplify(V*D*V.') S2 = simplify(inv(V*V.'))
S1 = [ 0, 0, 1, 2] [ 0, 1, 2, 0] [ 1, 2, 0, 0] [ 2, 0, 0, e] S2 = [ 0, 0, 0, 1] [ 0, 0, 1, 0] [ 0, 1, 0, 0] [ 1, 0, 0, 0]
Can carefully choosing the scaling of the eigenvectors be the basis for a sound numerical algorithm? I doubt it. We're still trying to compute something that is not unique, using factors that almost do not exist. It's pretty shaky.
- 범주:
- Matrices
댓글
댓글을 남기려면 링크 를 클릭하여 MathWorks 계정에 로그인하거나 계정을 새로 만드십시오.