/usr/share/doc/axiom-doc/hypertex/ocwmit18085lecture1.xhtml

<?xml version="1.0" encoding="UTF-8"?>
<html xmlns="http://www.w3.org/1999/xhtml" 
      xmlns:xlink="http://www.w3.org/1999/xlink"
      xmlns:m="http://www.w3.org/1998/Math/MathML">
 <head>
  <meta http-equiv="Content-Type" content="text/html" charset="us-ascii"/>
  <title>Axiom Documentation</title>
  <style>

   html {
     background-color: #ECEA81;
   }

   body { 
     margin: 0px;
     padding: 0px;
   }

   div.command { 
     color:red;
   }

   div.center {
     color:blue;
   }

   div.reset {
     visibility:hidden;
   }

   div.mathml { 
     color:blue;
   }

   input.subbut {
     background-color:#ECEA81;
     border: 0;
     color:green;
     font-family: "Courier New", Courier, monospace;
   }

   input.noresult {
     background-color:#ECEA81;
     border: 0;
     color:black;
     font-family: "Courier New", Courier, monospace;
   }

   span.cmd { 
     color:green;
     font-family: "Courier New", Courier, monospace;
   }

   pre {
     font-family: "Courier New", Courier, monospace;
   }
  </style>
 </head>
 <body>
  <div align="center"><img align="middle" src="doctitle.png"/></div>
  <hr/>
 Positive Definite Matrices K=A'CA
<hr/>
In applied mathematics we have 2 basic tasks:
<ul>
<li>Find the equations</li>
<li>Solve the equations</li>
</ul>
<h4>Positive Definite Matrices</h4>
Certain matrices occur frequently in applied math. These three
matrices (K,T,and M) are canonical examples.
We have 3 3x3 matrices, 
<pre>
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
</pre>
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:
<pre>
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))
</pre>
<pre>
t(n) == 
 M:=k(n)
 N:=M::Matrix(Fraction(Integer)) 
 qsetelt!(N,1,1,1) 
 N::SquareMatrix(n,Fraction(Integer))
</pre>
<pre>
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))
</pre>

K:=k(n) has a few key properties:
<ul>
<li> K is symmetric, that is K=K^T</li>
<li> K might be nonsingular, that is, it is invertible</li>
<li> K has a non-zero determinant</li>
<li> K is banded (main diagonal and neighbors)</li>
<li> K is tri-diagonal (main diagonal and nearest neighbors</li>
<li> K is extremely sparse</li>
<li> K has constant diagonals, (shift invariant, time invariant)</li>
<li> K is Toeplitz (constant diagonal, shows up in filters)</li>
<li> K is good for Fourier analysis</li>
</ul>

<h5>The inverse of T</h5>
If we look at the inverse of the T matrix we see:
<pre>
T^-1

        +3  2  1+
        |       |
        |2  2  1|
        |       |
        +1  1  1+
               Type: Matrix Fraction Integer
</pre>
Notice that these are all integers because the determinant of
this matrix is 1
<pre>
determinant T

     1
               Type: Fraction Integer

</pre>
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,
<pre>
[[-1,2,-1]]::MATRIX(INT)*[[3],[2],[1]]

     [0]
               Type: Matrix Integer

</pre>
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:
<pre>
 [[-1,2,-1]]::MATRIX(INT)*[[1],[1],[1]]

    [0]
               Type: Matrix Integer
</pre>
and the diagonal element of the unit matrix must be one. So
the second difference of the second column is:
<pre>
 [[-1,2,-1]]::MATRIX(INT)*[[2],[2],[1]]

    [1]
               Type: Matrix Integer
</pre>
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.

<br/>
<h5>The inverse of B</h5>
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:
<pre>
        + 1   - 1   0 + + x +    + 0 +
        |             | |   |    |   |
        |- 1   2   - 1| | y | =  | 0 |
        |             | |   |    |   |
        + 0   - 1   1 + + z +    + 0 +

</pre>
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
<pre>
   B x = 0
</pre>
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
<pre>
        + 1   - 1   0 + + x +    + 0 +
        |             | |   |    |   |
        |- 1   2   - 1| | y | =  | 0 |
        |             | |   |    |   |
        + 0   - 1   1 + + z +    + 0 +

</pre>
we can think about this as 
<pre>
x*(first column) + y*(second column) + z*(third column).
</pre>
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 
<pre>
 (first row)*(vector)
 (second row)*(vector)
 (third row)*(vector)
</pre>

<br/>
<h5>The inverse of K</h5>
Now we consider the K matrix we see the inverse
<pre>
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)
</pre>
We can take the determinant of k 
<pre>
determinant K

    4
               Type: Fraction Integer
</pre>
Thus there is a constant 1/4 which can be factored out
<pre>
4*kinv

         +3  2  1+
         |       |
         |2  4  2|
         |       |
         +1  2  3+
               Type: SquareMatrix(3,Fraction Integer)
</pre>
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
<pre>
 K u = f
</pre>
by elimination which implies multiplying and subtracting rows.
<pre>
       K    u  =  f    ==>   U     u  =    f
</pre>                                        
For the 2x2 case we see:
<pre>
                             +2  -1+        +  f1  +
    +2  -1+  +x+   +f1+      |     |  +x+   |      |
    |     |  | | = |  |  ==> |    3|  | | = |   1  |
    +-1  2+  +y+   +f2+      |0   -|  +y+   |f2+-f1|
                             +    2+        +   2  +


</pre>
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.

<hr/>
<h5>Generalizing the matrix pivot operations</h5>
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
<pre>
      E     K    =   U
</pre>
In detail we see
<pre>

    +1  0+            +2  -1+
    |    |  +2  -1+   |     |
    |1   |  |     | = |    3|
    |-  1|  +-1  2+   |0   -|
    +2   +            +    2+

</pre>
We wish to rewrite this as
<pre>
       K = L U 
</pre>

<hr/>
<h5>The big 4 solve operations in Linear Algebra</h5>
<ol>
<li>Elimination</li>
<li>Gram-Schmidt Orthoginalization</li>
<li>Eigenvalues</li>
<li>Singular Value Decomposition</li>
</ol>
Each of these operations is described by a factorization of K.
Elimination is written 
<pre>
  K = L U
</pre>
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
<pre>

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
</pre>
Given the matrix operations above we had
<pre>
      E       K   =   U

    +1  0+            +2  -1+
    |    |  +2  -1+   |     |
    |1   |  |     | = |    3|
    |-  1|  +-1  2+   |0   -|
    +2   +            +    2+

</pre>
and the inverse of E is the same matrix with a minus sign in
the second row, thus:
<pre>
        +  1  0+ 
   -1   |      | 
  E   = |  1   | = L 
        |- -  1| 
        +  2   + 

</pre>

<hr/>
<h5>Making the matrices symmetric</h5>
We would like to preserve the symmetry property which we can
do with a further decomposition of LU as follows:
<pre>
      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+

</pre>
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
<pre>
           T
    L  D  L
</pre>
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
<pre>
            T  T     TT  T   T        T T        T
  (  L  D  L  )   = L   D   L   =  L D L  = L D L
</pre>
<hr/>
<h5>Positive Definite Matrices</h5>
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
<pre>
   positive definite      == all pivots >  0
   positive semi-definite == all pivots >= 0
</pre>
When all the pivots are positive then all the eigenvalues are positive.

So a positive definite matrix K and any non-zero vector X
<pre>
    T
   X  K X  > 0
</pre>
X transpose is just a row and X is just a column.

 </body>
</html>
axiom-hypertex-data 20120501-8 / usr / share / doc / axiom-doc / hypertex / ocwmit18085lecture1.xhtml
