Difference between revisions of "Software tutorial/Matrix operations"

From Process Model Formulation and Solution: 3E4
Jump to navigation Jump to search
m (Created page with "== Matrix operations == In this section we will see how to calculate the LU decomposition, matrix inverse, norms and condition number of a matrix. Let A = \(\left(\begin{arra...")
 
m
 
(7 intermediate revisions by the same user not shown)
Line 1: Line 1:
== Matrix operations ==
{{Navigation|Book=Software tutorial|previous=Vectors and arrays|current=Tutorial index|next=Functions as objects}}
__TOC__
== LU decomposition ==


In this section we will see how to calculate the LU decomposition, matrix inverse, norms and condition number of a matrix. 
Let \(A = \left(\begin{array}{cccc}     4 & 0  & 1 \\ 0  & 6 & 2 \\ 3.2  & 7 & 4 \end{array} \right)\) and let \(b = \left(\begin{array}{cccc} 9.1 \\ 2  \\ 0 \end{array} \right) \)
 
Let A = \(\left(\begin{array}{cccc} & 0  & 0 \\ 0  & 1 & 0 \\ 0 & 0 & 1 \end{array} \right)\) and let B = \(\left(\begin{array}{cccc}     4 & 0  & 1 \\ 0  & 6  & 2  \\ 3.2  & 7  &  4 \end{array} \right)\)


{| class="wikitable"
{| class="wikitable"
Line 12: Line 12:
| width="50%" valign="top" |
| width="50%" valign="top" |
<syntaxhighlight lang="matlab">
<syntaxhighlight lang="matlab">
A = eye(3);
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];
B = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];


% LU decomposition of B
% LU decomposition of A
[L,U,P] = lu(B)
[L,U,P] = lu(A)
L =
L =
     1.0000        0        0
     1.0000        0        0
Line 30: Line 29:
         0        1        0
         0        1        0
</syntaxhighlight>
</syntaxhighlight>
The <tt>P</tt> matrix shows that rows 2 and rows 3 were interchanged during the Gauss elimination pivoting steps.
| width="50%" valign="top" |
| width="50%" valign="top" |
<syntaxhighlight lang="matlab">
<syntaxhighlight lang="python">
import numpy as np
import numpy as np
from numpy.linalg import *
from scipy.linalg import *
from scipy.linalg import *


A = np.eye(3)
A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])
B = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])
LU, P = lu_factor(A)
LU, P = lu_factor(B)
 
LU
>>> LU
array([[ 4.        ,  0.        ,  1.        ],
array([[ 4.        ,  0.        ,  1.        ],
       [ 0.8      ,  7.        ,  3.2      ],
       [ 0.8      ,  7.        ,  3.2      ],
       [ 0.        ,  0.85714286, -0.74285714]])
       [ 0.        ,  0.85714286, -0.74285714]])


P
>>> P
array([0, 2, 2])
array([0, 2, 2])
</syntaxhighlight>
</syntaxhighlight>


The LU matrices are overlayed in the output.  Notice that the results are the same as MATLAB's, except we do not need to store the diagonal set of ones from the L matrix.
The P vector is also a bit cryptic - it should be read from left to right, in order.  The way to read it is that row \(i\) of matrix \(A\) was interchanged with row <tt>P[i]</tt> in \(A\), where <tt>i</tt> is the index into vector <tt>P</tt>.
* row 0 in \(A\) was interchanged with row 0: i.e. no pivoting for the first row
* row 1 in \(A\) was interchanged with row 2: i.e. rows 2 and 3 in \(A\) were swapped
* row 2 in \(A\) was interchanged with row 2: i.e. no further pivoting.
|}
Once we have \(L, U\) and \(P\) we can solve the system of equations for a given vector, \(b\).  It takes a bit more work to do this in MATLAB than with Python.  Recall from class that \(A = LU\), but with a permutation or pivot matrix, it becomes \(PA = LU\), or alternatively: \(A = P^{-1}LU\)
\[ \begin{align}  Ax &= b \\ P^{-1}LUx &= b \\ \text{Let}\qquad Ux &= y \\ \text{Then solve}\qquad \left(P^{-1}L\right)y &= b \qquad \text{for}\,\, y \\ \text{Finally, solve}\qquad Ux &= y \qquad \text{for}\,\, x \end{align}  \]
{| class="wikitable"
|-
! MATLAB
! Python
|-
| width="50%" valign="top" |
<syntaxhighlight lang="matlab">
b = [9.1, 2, 0]';
[L,U,P] = lu(A);
y = mldivide(inv(P)*L, b);
x = mldivide(U, y);
>> x
x =
    5.0481
    4.0308
  -11.0923
% Verify the solution:
>> A*x - b
ans =
  1.0e-14 *
    0.1776
        0
        0
</syntaxhighlight>
The small values of the solution error indicate that our LU solution is accurate.  You should compare the LU solution with solution from <tt>mldivide(A, b)</tt>.
| width="50%" valign="top" |
<syntaxhighlight lang="python">
LU, P = lu_factor(A)
b = np.array([9.1, 2, 0])
x = lu_solve((LU,P), b)
>>> x
array([  5.04807692,  4.03076923, -11.09230769])
# Verify the solution:
>>> np.dot(A, x) - b
array([  1.77635684e-15,  0.00000000e+00,  0.00000000e+00])
</syntaxhighlight>
|}
== Matrix determinant ==
One way to test if a matrix is invertable is to calculate the determinant.  If it is non-zero, then \(A\) is invertable.
Let \(A = \left(\begin{array}{cccc}    4 & 0  & 1 \\ 0  & 6  & 2  \\ 3.2  & 7  &  4 \end{array} \right)\)
{| class="wikitable"
|-
! MATLAB
! Python
|-
| width="50%" valign="top" |
<syntaxhighlight lang="matlab">
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];
>> det(A)
ans =
  20.8000
</syntaxhighlight>
| width="50%" valign="top" |
<syntaxhighlight lang="python">
import numpy as np
from numpy.linalg import *
from scipy.linalg import *
A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])
>>> det(A)
20.800000000000011
</syntaxhighlight>
|}
== Matrix inverse ==
The matrix inverse is defined so that \(A A^{-1} = I\).  Let \(A = \left(\begin{array}{cccc}    4 & 0  & 1 \\ 0  & 6  & 2  \\ 3.2  & 7  &  4 \end{array} \right)\)
{| class="wikitable"
|-
! MATLAB
! Python
|-
| width="50%" valign="top" |
<syntaxhighlight lang="matlab">
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];
>> inv(A)
ans =
    0.4808    0.3365  -0.2885
    0.3077    0.6154  -0.3846
  -0.9231  -1.3462    1.1538
>> A * inv(A)
ans =
    1.0000        0        0
  -0.0000    1.0000    0.0000
        0        0    1.0000
</syntaxhighlight>
| width="50%" valign="top" |
<syntaxhighlight lang="python">
import numpy as np
from numpy.linalg import *
from scipy.linalg import *
A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])
>>> inv(A)
array([[ 0.48076923,  0.33653846, -0.28846154],
      [ 0.30769231,  0.61538462, -0.38461538],
      [-0.92307692, -1.34615385,  1.15384615]])
>>> np.dot(A, inv(A))
array([[  1.00000000e+00,  0.00000000e+00,  0.00000000e+00],
      [ -2.22044605e-16,  1.00000000e+00,  4.44089210e-16],
      [  0.00000000e+00,  0.00000000e+00,  1.00000000e+00]])
</syntaxhighlight>
Notice the small, essentially zero elements in the off-diagonal elements.
|}
== Matrix and vector norms ==
{| class="wikitable"
|-
! MATLAB
! Python
|-
| width="50%" valign="top" |
<syntaxhighlight lang="matlab">
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];
format long
% The default norm in MATLAB = max(svd(A))
>> norm(A)
ans =
  10.627892098461123
>> max(svd(A))
ans =
  10.627892098461123
% Other norms available:
>> norm(A, 'fro')  % Frobenius norm: square root of
                  % the sum of squares of all entries in A
ans =
  11.499565209172042
% Other norms
% norm(A, 1)  : column sum norm: 13.0
% norm(A, inf) : row sum norm: 14.199999999999999
</syntaxhighlight>
| width="50%" valign="top" |
<syntaxhighlight lang="python">
import numpy as np
from numpy.linalg import *
from scipy.linalg import *
A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])
# The default norm in Python = Frobenius norm
>>> norm(A)
ans =
  11.499565209172042
>>> norm(A, 'fro')
ans =
  11.499565209172042
# Other norms available:
>>> norm(A, 2)  % 2-norm: maximum singular value
ans =
  10.627892098461123
# Other norms
# norm(A, 1)      : column sum norm: 13.0
# norm(A, np.inf) : row sum norm: 14.199999999999999
</syntaxhighlight>
|}
== Matrix conditioning ==
The condition number was defined as \({\rm cond}(A) \stackrel{\Delta}{=} \|A\|\times\|A^{-1}\|\).  Any norm may be used, but if left unstated, then you should assume the 2-norm was used (maximum singular value).
{| class="wikitable"
|-
! MATLAB
! Python
|-
| width="50%" valign="top" |
<syntaxhighlight lang="matlab">
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];
format long
>> cond(A)
ans =
  23.724794742061690
% Condition numbers, calculated using other norms are also available
>> help cond
COND  Condition number with respect to inversion.
    COND(X) returns the 2-norm condition number (the ratio of the
    largest singular value of X to the smallest).  Large condition
    numbers indicate a nearly singular matrix.
    COND(X,P) returns the condition number of X in P-norm:
      NORM(X,P) * NORM(INV(X),P).
    where P = 1, 2, inf, or 'fro'.
</syntaxhighlight>
| width="50%" valign="top" |
<syntaxhighlight lang="python">
import numpy as np
from numpy.linalg import *
from scipy.linalg import *
A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])
>>> cond(A)
23.72479474206169
# Condition numbers, calculated using other norms are also available
>>> help(cond)
cond(x, p=None)
    Compute the condition number of a matrix.
   
    This function is capable of returning the condition number using
    one of seven different norms, depending on the value of `p` (see
    Parameters below).
   
    Parameters
    ----------
    x : array_like, shape (M, N)
        The matrix whose condition number is sought.
    p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional
        Order of the norm:
   
        =====  ============================
        p      norm for matrices
        =====  ============================
        None  2-norm, computed directly using the ``SVD``
        'fro'  Frobenius norm
        inf    max(sum(abs(x), axis=1))
        -inf  min(sum(abs(x), axis=1))
        1      max(sum(abs(x), axis=0))
        -1    min(sum(abs(x), axis=0))
        2      2-norm (largest sing. value)
        -2    smallest singular value
        =====  ============================
</syntaxhighlight>
|}
|}
{{Navigation|Book=Software tutorial|previous=Vectors and arrays|current=Tutorial index|next=Functions as objects}}

Latest revision as of 16:05, 18 October 2010

← Vectors and arrays (previous step) Tutorial index Next step: Functions as objects →

LU decomposition

Let \(A = \left(\begin{array}{cccc} 4 & 0 & 1 \\ 0 & 6 & 2 \\ 3.2 & 7 & 4 \end{array} \right)\) and let \(b = \left(\begin{array}{cccc} 9.1 \\ 2 \\ 0 \end{array} \right) \)

MATLAB Python
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];

% LU decomposition of A
[L,U,P] = lu(A)
L =
    1.0000         0         0
    0.8000    1.0000         0
         0    0.8571    1.0000
U =
    4.0000         0    1.0000
         0    7.0000    3.2000
         0         0   -0.7429
P =
         1         0         0
         0         0         1
         0         1         0

The P matrix shows that rows 2 and rows 3 were interchanged during the Gauss elimination pivoting steps.

import numpy as np
from numpy.linalg import *
from scipy.linalg import *

A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])
LU, P = lu_factor(A)

>>> LU
array([[ 4.        ,  0.        ,  1.        ],
       [ 0.8       ,  7.        ,  3.2       ],
       [ 0.        ,  0.85714286, -0.74285714]])

>>> P
array([0, 2, 2])

The LU matrices are overlayed in the output. Notice that the results are the same as MATLAB's, except we do not need to store the diagonal set of ones from the L matrix.

The P vector is also a bit cryptic - it should be read from left to right, in order. The way to read it is that row \(i\) of matrix \(A\) was interchanged with row P[i] in \(A\), where i is the index into vector P.

  • row 0 in \(A\) was interchanged with row 0: i.e. no pivoting for the first row
  • row 1 in \(A\) was interchanged with row 2: i.e. rows 2 and 3 in \(A\) were swapped
  • row 2 in \(A\) was interchanged with row 2: i.e. no further pivoting.

Once we have \(L, U\) and \(P\) we can solve the system of equations for a given vector, \(b\). It takes a bit more work to do this in MATLAB than with Python. Recall from class that \(A = LU\), but with a permutation or pivot matrix, it becomes \(PA = LU\), or alternatively: \(A = P^{-1}LU\)

\[ \begin{align} Ax &= b \\ P^{-1}LUx &= b \\ \text{Let}\qquad Ux &= y \\ \text{Then solve}\qquad \left(P^{-1}L\right)y &= b \qquad \text{for}\,\, y \\ \text{Finally, solve}\qquad Ux &= y \qquad \text{for}\,\, x \end{align} \]

MATLAB Python
b = [9.1, 2, 0]';
[L,U,P] = lu(A);
y = mldivide(inv(P)*L, b);
x = mldivide(U, y);
>> x
x =
    5.0481
    4.0308
  -11.0923

% Verify the solution:
>> A*x - b
ans =
   1.0e-14 *

    0.1776
         0
         0

The small values of the solution error indicate that our LU solution is accurate. You should compare the LU solution with solution from mldivide(A, b).

LU, P = lu_factor(A)
b = np.array([9.1, 2, 0])
x = lu_solve((LU,P), b)

>>> x
array([  5.04807692,   4.03076923, -11.09230769])

# Verify the solution:
>>> np.dot(A, x) - b
array([  1.77635684e-15,   0.00000000e+00,   0.00000000e+00])

Matrix determinant

One way to test if a matrix is invertable is to calculate the determinant. If it is non-zero, then \(A\) is invertable.

Let \(A = \left(\begin{array}{cccc} 4 & 0 & 1 \\ 0 & 6 & 2 \\ 3.2 & 7 & 4 \end{array} \right)\)

MATLAB Python
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];
>> det(A)
ans =
   20.8000
import numpy as np
from numpy.linalg import *
from scipy.linalg import *

A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])
>>> det(A)
20.800000000000011

Matrix inverse

The matrix inverse is defined so that \(A A^{-1} = I\). Let \(A = \left(\begin{array}{cccc} 4 & 0 & 1 \\ 0 & 6 & 2 \\ 3.2 & 7 & 4 \end{array} \right)\)

MATLAB Python
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];
>> inv(A)
ans =
    0.4808    0.3365   -0.2885
    0.3077    0.6154   -0.3846
   -0.9231   -1.3462    1.1538

>> A * inv(A)
ans =
    1.0000         0         0
   -0.0000    1.0000    0.0000
         0         0    1.0000
import numpy as np
from numpy.linalg import *
from scipy.linalg import *

A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])
>>> inv(A)
array([[ 0.48076923,  0.33653846, -0.28846154],
       [ 0.30769231,  0.61538462, -0.38461538],
       [-0.92307692, -1.34615385,  1.15384615]])

>>> np.dot(A, inv(A))
array([[  1.00000000e+00,   0.00000000e+00,   0.00000000e+00],
       [ -2.22044605e-16,   1.00000000e+00,   4.44089210e-16],
       [  0.00000000e+00,   0.00000000e+00,   1.00000000e+00]])

Notice the small, essentially zero elements in the off-diagonal elements.

Matrix and vector norms

MATLAB Python
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];
format long


% The default norm in MATLAB = max(svd(A))
>> norm(A)
ans =
  10.627892098461123
>> max(svd(A))
ans =
  10.627892098461123

% Other norms available: 
>> norm(A, 'fro')  % Frobenius norm: square root of
                   % the sum of squares of all entries in A
ans =
  11.499565209172042

% Other norms
% norm(A, 1)   : column sum norm: 13.0
% norm(A, inf) : row sum norm: 14.199999999999999
import numpy as np
from numpy.linalg import *
from scipy.linalg import *

A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])

# The default norm in Python = Frobenius norm
>>> norm(A)
ans =
  11.499565209172042
>>> norm(A, 'fro')
ans =
  11.499565209172042

# Other norms available: 
>>> norm(A, 2)  % 2-norm: maximum singular value
ans =
  10.627892098461123

# Other norms
# norm(A, 1)      : column sum norm: 13.0
# norm(A, np.inf) : row sum norm: 14.199999999999999


Matrix conditioning

The condition number was defined as \({\rm cond}(A) \stackrel{\Delta}{=} \|A\|\times\|A^{-1}\|\). Any norm may be used, but if left unstated, then you should assume the 2-norm was used (maximum singular value).

MATLAB Python
A = [4, 0, 1; 0, 6, 2; 3.2, 7, 4];
format long

>> cond(A)
ans =
  23.724794742061690



% Condition numbers, calculated using other norms are also available

>> help cond
 COND   Condition number with respect to inversion.
    COND(X) returns the 2-norm condition number (the ratio of the
    largest singular value of X to the smallest).  Large condition
    numbers indicate a nearly singular matrix.
 
    COND(X,P) returns the condition number of X in P-norm:
 
       NORM(X,P) * NORM(INV(X),P). 
 
    where P = 1, 2, inf, or 'fro'.
import numpy as np
from numpy.linalg import *
from scipy.linalg import *

A = np.array([[4, 0, 1],[0, 6, 2],[3.2, 7, 4]])

>>> cond(A)
23.72479474206169

# Condition numbers, calculated using other norms are also available

>>> help(cond)
cond(x, p=None)
    Compute the condition number of a matrix.
    
    This function is capable of returning the condition number using
    one of seven different norms, depending on the value of `p` (see
    Parameters below).
    
    Parameters
    ----------
    x : array_like, shape (M, N)
        The matrix whose condition number is sought.
    p : {None, 1, -1, 2, -2, inf, -inf, 'fro'}, optional
        Order of the norm:
    
        =====  ============================
        p      norm for matrices
        =====  ============================
        None   2-norm, computed directly using the ``SVD``
        'fro'  Frobenius norm
        inf    max(sum(abs(x), axis=1))
        -inf   min(sum(abs(x), axis=1))
        1      max(sum(abs(x), axis=0))
        -1     min(sum(abs(x), axis=0))
        2      2-norm (largest sing. value)
        -2     smallest singular value
        =====  ============================
← Vectors and arrays (previous step) Tutorial index Next step: Functions as objects →