Positive Definite Matrices K=A'CA
In applied mathematics we have 2 basic tasks:
- Find the equations
- Solve the equations
Positive Definite Matrices
Certain matrices occur frequently in applied math. These three
matrices (K,T,and M) are canonical examples.
We have 3 3x3 matrices,
K:Matrix(Integer):=[[2,-1,0],[-1,2,-1],[0,-1,2]]
+ 2 - 1 0 +
| |
|- 1 2 - 1|
| |
+ 0 - 1 2 +
Type: Matrix Integer
T:Matrix(Integer):=[[1,-1,0],[-1,2,-1],[0,-1,2]]
+ 1 - 1 0 +
| |
|- 1 2 - 1|
| |
+ 0 - 1 2 +
Type: Matrix Integer
B:Matrix(Integer):=[[1,-1,0],[-1,2,-1],[0,-1,1]]
+ 1 - 1 0 +
| |
|- 1 2 - 1|
| |
+ 0 - 1 1 +
Type: Matrix Integer
These matrices are similar and can be generalized to square matrices
of order N, with n x n elements. All of these matrices have the same
element along the diagonal. T (aka Top) differs from K in the first row.
B (aka Both) differs from K in the first and last row. These represent
different boundary conditions in the problem.
We can create K(n), T(n) and B(n) with the following commands:
k(n) ==
M := diagonalMatrix([2 for i in 1..n])
for i in 1..n-1 repeat M(i,i+1):=-1
for i in 1..n-1 repeat M(i+1,i):=-1
M::SquareMatrix(n,Fraction(Integer))
t(n) ==
M:=k(n)
N:=M::Matrix(Fraction(Integer))
qsetelt!(N,1,1,1)
N::SquareMatrix(n,Fraction(Integer))
b(n) ==
M:=k(n)
N:=M::Matrix(Fraction(Integer))
qsetelt!(N,1,1,1)
qsetelt!(N,n,n,1)
N::SquareMatrix(n,Fraction(Integer))
K:=k(n) has a few key properties:
- K is symmetric, that is K=K^T
- K might be nonsingular, that is, it is invertible
- K has a non-zero determinant
- K is banded (main diagonal and neighbors)
- K is tri-diagonal (main diagonal and nearest neighbors
- K is extremely sparse
- K has constant diagonals, (shift invariant, time invariant)
- K is Toeplitz (constant diagonal, shows up in filters)
- K is good for Fourier analysis
The inverse of T
If we look at the inverse of the T matrix we see:
T^-1
+3 2 1+
| |
|2 2 1|
| |
+1 1 1+
Type: Matrix Fraction Integer
Notice that these are all integers because the determinant of
this matrix is 1
determinant T
1
Type: Fraction Integer
We can check that this matrix is the inverse of T.
When computing the inverse the row pattern [-1 2 -1] is a
``second difference''. The first column of the inverse matrix
is [3 2 1] which is linear. When we take the second difference
of a linear object we should get 0. Thus,
[[-1,2,-1]]::MATRIX(INT)*[[3],[2],[1]]
[0]
Type: Matrix Integer
The third column of the T matrix is linear and constant. If we
take the second difference of that we also find it is zero:
[[-1,2,-1]]::MATRIX(INT)*[[1],[1],[1]]
[0]
Type: Matrix Integer
and the diagonal element of the unit matrix must be one. So
the second difference of the second column is:
[[-1,2,-1]]::MATRIX(INT)*[[2],[2],[1]]
[1]
Type: Matrix Integer
So these simple checks show that we're getting the correct
row and column values for the identity matrix by multiplying
T times its inverse.
The inverse of B
If we look for the inverse of the B matrix we can observe
that the rows sum to zero which implies that it is not
invertible. Thus it is singular.
K and T are positive definite. B is only positive semi-definite.
If we can find a vector that it takes to zero, that is if we can
solve for x,y,z in:
+ 1 - 1 0 + + x + + 0 +
| | | | | |
|- 1 2 - 1| | y | = | 0 |
| | | | | |
+ 0 - 1 1 + + z + + 0 +
The constant vector [1 1 1] solves this equation. When
the rows sum to zero we are adding each row by a constant
and thus we add each row times the constant one and we
get zeros. If the matrix takes some vector to zero it
cannot have an inverse since if
B x = 0
and x is not zero. If B had an inverse only x=0 would
solve the equation. Since x=1 solves the equation B has
no inverse. The vector x is in the nullspace of B. In
fact any constant vector, e.g. [3 3 3] is in the nullspace.
Thus the nullspace of B is cx for any constant c.
When doing matrix multiplication one way to think about the
work is to consider the problem by columns. Thus in the
multiplication
+ 1 - 1 0 + + x + + 0 +
| | | | | |
|- 1 2 - 1| | y | = | 0 |
| | | | | |
+ 0 - 1 1 + + z + + 0 +
we can think about this as
x*(first column) + y*(second column) + z*(third column).
and for the constant vector [1 1 1] this means that we
just need to sum the columns.
Alternatively this can be computed by thinking of the
multiplication as
(first row)*(vector)
(second row)*(vector)
(third row)*(vector)
The inverse of K
Now we consider the K matrix we see the inverse
K
+ 2 - 1 0 +
| |
|- 1 2 - 1|
| |
+ 0 - 1 2 +
Type: SquareMatrix(3,Fraction Integer)
kinv:=K^-1
+3 1 1+
|- - -|
|4 2 4|
| |
|1 1|
|- 1 -|
|2 2|
| |
|1 1 3|
|- - -|
+4 2 4+
Type: SquareMatrix(3,Fraction Integer)
We can take the determinant of k
determinant K
4
Type: Fraction Integer
Thus there is a constant 1/4 which can be factored out
4*kinv
+3 2 1+
| |
|2 4 2|
| |
+1 2 3+
Type: SquareMatrix(3,Fraction Integer)
Notice that the inverse is a symmetric matrix but not tri-diagonal.
The inverse is not a sparse matrix so much more computation would
be involved when using the inverse.
In order to solve the system
K u = f
by elimination which implies multiplying and subtracting rows.
K u = f ==> U u = f
For the 2x2 case we see:
+2 -1+ + f1 +
+2 -1+ +x+ +f1+ | | +x+ | |
| | | | = | | ==> | 3| | | = | 1 |
+-1 2+ +y+ +f2+ |0 -| +y+ |f2+-f1|
+ 2+ + 2 +
By multiplying row1 by 1/2 and adding it to row2 we create an
upper triangular matrix U. Since we chose K(1,1), the number 2
is called the first pivot. K(2,2), the number 3/2, is called
the second pivot.
For K 2x2 above is symmetric and invertible (since the pivots
are all non-zero).
For the K 3x3 case the pivots are 2, 3/2, and 4/3. (The next pivots
would be 5/4, 6/5, etc. for larger matrices).
For the T 3x3 case the pivots are 1, 1, and 1.
For the B 3x3 case the third pivot would be zero.
Generalizing the matrix pivot operations
For the 2x2 case we see contruct an elimination matrix E which we can use
to pre-multipy by K to give us the upper triangular matrix U
E K = U
In detail we see
+1 0+ +2 -1+
| | +2 -1+ | |
|1 | | | = | 3|
|- 1| +-1 2+ |0 -|
+2 + + 2+
We wish to rewrite this as
K = L U
The big 4 solve operations in Linear Algebra
- Elimination
- Gram-Schmidt Orthoginalization
- Eigenvalues
- Singular Value Decomposition
Each of these operations is described by a factorization of K.
Elimination is written
K = L U
where L is lower triangular and U is upper triangular.
Thus we need a matrix L which when multiplied by U gives K.
The required matrix is the inverse of the E matrix above since
1) E K = U
-1 -1
2) E E K = E U
-1
3) I K = E U
-1
4) but L = E
5) so K = L U
Given the matrix operations above we had
E K = U
+1 0+ +2 -1+
| | +2 -1+ | |
|1 | | | = | 3|
|- 1| +-1 2+ |0 -|
+2 + + 2+
and the inverse of E is the same matrix with a minus sign in
the second row, thus:
+ 1 0+
-1 | |
E = | 1 | = L
|- - 1|
+ 2 +
Making the matrices symmetric
We would like to preserve the symmetry property which we can
do with a further decomposition of LU as follows:
L U = L D U'
+ 1 0+ +2 -1+ + 1 0+ +2 0+ +1 1+
| | | | | | | | | - -|
| 1 | | 3| = | 1 | | 3| | 2|
|- - 1| |0 -| |- - 1| |0 -| | |
+ 2 + + 2+ + 2 + + 2+ +0 1+
So now we have 3 matrices; L is the lower triangular,
D is symmetric and contains the pivots, and U' is upper triangular and
is the transpose of the lower. So the real form we have is
T
L D L
This result will always be symmetric. We can check this by taking
its transpose. If we get the same matrix we must have a symmetric
matrix. So the transpose of
T T TT T T T T T
( L D L ) = L D L = L D L = L D L
Positive Definite Matrices
There are several ways to recognize a positive definite matrix.
First, it must be symmetric. The "positive" aspect comes from
the pivots, all of which must be positive. Note that T is also
positive definite. B is positive semi-definite because one of
the pivots is zero. So
positive definite == all pivots > 0
positive semi-definite == all pivots >= 0
When all the pivots are positive then all the eigenvalues are positive.
So a positive definite matrix K and any non-zero vector X
T
X K X > 0
X transpose is just a row and X is just a column.
