MAT 232 Chapter Five

506 days ago by kcrisman

Similarity and Diagonalizability

Here is an example of a diagonalizable matrix and what happens.  Can you "see" how it does spread out evenly in precisely two directions?

t = var('t') @interact def _(A=matrix(RDF,[[4,-2],[1,1]]),auto_update=False): ID = matrix(RDF,[[1,0],[0,1]]) def makeface(M): pll=M*vector((-0.5,0.5)) plr=M*vector((-0.3,0.5)) prl=M*vector((0.3,0.5)) prr=M*vector((0.5,0.5)) left_eye=line([pll,plr])+point(pll,size=5)+point(plr,size=5) right_eye=line([prl,prr],color='green')+point(prl,size=5,color='green')+point(prr,size=5,color='green') mouth=parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t, -0.5, 0),color='red')+parametric_plot(M*vector([t, -0.15*sin(2*pi*t)-0.5]), (t,0,0.5),color='orange') face=parametric_plot(M*vector([cos(t),sin(t)]), (t,0,pi/2),color='black')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi/2,pi),color='lavender')+parametric_plot(M*vector([cos(t),sin(t)]), (t,pi,3*pi/2),color='cyan')+parametric_plot(M*vector([cos(t),sin(t)]),(t,3*pi/2,2*pi),color='sienna') return right_eye+left_eye+face+mouth pretty_print(html('smiley guy, then transformed by $A$')) G = makeface(ID)+makeface(A) G.show(aspect_ratio=1,figsize=5) 
       

Click to the left again to hide and once more to show the dynamic interactive window

 
       

Powers of diagonalizable matrices!  You have to get a pretty high power, but eventually diagonalization is faster.

M = matrix(3,[0,1,1,1,0,1,1,1,0]) 
       
timeit('M^100000') 
       
25 loops, best of 3: 9.74 ms per loop
25 loops, best of 3: 9.74 ms per loop
D,P=M.jordan_form(transformation=True); Pinv=P^-1; D = D.change_ring(ZZ) 
       
timeit('D,P=M.jordan_form(transformation=True); Pinv=P^-1') 
       
25 loops, best of 3: 39.3 ms per loop
25 loops, best of 3: 39.3 ms per loop
timeit('Pinv*D^100000*P') 
       
125 loops, best of 3: 1.27 ms per loop
125 loops, best of 3: 1.27 ms per loop
 
       
 
       
 
       

Recall we did some data from orthogonalization and least squares.

 
       

Orthogonal diagonalization?

M = matrix(QQ,[[1,1,3],[1,3,1],[3,1,1]]); show(M) 
       
M.charpoly() 
       
x^3 - 5*x^2 - 4*x + 20
x^3 - 5*x^2 - 4*x + 20
M.eigenvalues() 
       
[5, 2, -2]
[5, 2, -2]
M.is_diagonalizable() 
       
True
True
M.jordan_form(transformation=True) 
       
(
[ 5| 0| 0]            
[--+--+--]            
[ 0| 2| 0]  [ 1  1  1]
[--+--+--]  [ 1 -2  0]
[ 0| 0|-2], [ 1  1 -1]
)
(
[ 5| 0| 0]            
[--+--+--]            
[ 0| 2| 0]  [ 1  1  1]
[--+--+--]  [ 1 -2  0]
[ 0| 0|-2], [ 1  1 -1]
)
P = matrix([[1,1,1],[1,-2,0],[1,1,-1]]) 
       
P^-1*M*P 
       
[ 5  0  0]
[ 0  2  0]
[ 0  0 -2]
[ 5  0  0]
[ 0  2  0]
[ 0  0 -2]

Note that the columns of the change of basis matrix form an orthogonal basis for $\mathbb{R}^3$!

 
       

One application of orthogonal diagonalization we have already seen, or at least a close cousin:

import pylab A_image = pylab.mean(pylab.imread(DATA + 'Jenks.png'), 2) @interact def svd_image(i=(5,(1..100))): u,s,v = pylab.linalg.svd(A_image) A = sum(s[j]*pylab.outer(u[0:,j],v[j,0:]) for j in range(i)) g = graphics_array([[matrix_plot(A),matrix_plot(A_image)]]) show(g,axes=False, figsize=[8,10]) pretty_print(html('<h2>Compressed using %s eigenvalues</h2>'%i)) 
       

Click to the left again to hide and once more to show the dynamic interactive window

 
       

But there are other very important applications too.

Let's use some data about liver problems from http://archive.ics.uci.edu/ml/machine-learning-databases/liver-disorders/bupa.data

1. Title: BUPA liver disorders

2. Source information:
   -- Creators: BUPA Medical Research Ltd.
   -- Donor: Richard S. Forsyth
             8 Grosvenor Avenue
             Mapperley Park
             Nottingham NG3 5DX
             0602-621676
   -- Date: 5/15/1990

3. Past usage: 
   -- None known other than what is shown in the PC/BEAGLE User's Guide
      (written by Richard S. Forsyth).

4. Relevant information:
   -- The first 5 variables are all blood tests which are thought
      to be sensitive to liver disorders that might arise from
      excessive alcohol consumption.  Each line in the bupa.data file
      constitutes the record of a single male individual.
   -- It appears that drinks>5 is some sort of a selector on this database.
      See the PC/BEAGLE User's Guide for more information.

5. Number of instances: 345

6. Number of attributes: 7 overall

7. Attribute information:
   1. mcv	mean corpuscular volume
   2. alkphos	alkaline phosphotase
   3. sgpt	alamine aminotransferase
   4. sgot 	aspartate aminotransferase
   5. gammagt	gamma-glutamyl transpeptidase
   6. drinks	number of half-pint equivalents of alcoholic beverages
                drunk per day
   7. selector  field used to split data into two sets

8. Missing values: none
import csv csvfile = open(DATA+'liver.csv','rb') R = csv.reader(csvfile, delimiter=',') data = list(R) csvfile.close() 
       
len(data) 
       
345
345
data = [[QQ(RDF(x)) for x in row][:-1] for row in data] 
       
M = matrix(data); M 
       
345 x 6 dense matrix over Rational Field (use the '.str()' method to see
the entries)
345 x 6 dense matrix over Rational Field (use the '.str()' method to see the entries)
M = M.transpose() 
       
m = vector([mean(row) for row in M]); m 
       
(6221/69, 1607/23, 2098/69, 2834/115, 13208/345, 1192/345)
(6221/69, 1607/23, 2098/69, 2834/115, 13208/345, 1192/345)
B = matrix(RDF,[vector(d)-m for d in M.columns()]).transpose() 
       
S = 1/(B.ncols())*B*B.transpose() 
       
       
[  19.7282083595883  3.588909892879645 12.781684520058807 
8.381474480151232  38.70544003360641  4.628901491283342]
[ 3.588909892879647  335.6612476370509 27.203654694391926
26.892627599243863  95.61386263390044  6.155009451795841]
[12.781684520058807 27.203654694391926  379.6266330602812 
144.8374291115312 384.48762864944314 13.432724217601342]
[ 8.381474480151232 26.892627599243863  144.8374291115312
101.00042848141132 207.84909892879648  9.365141776937618]
[ 38.70544003360641  95.61386263390044 384.48762864944314
207.84909892879648 1536.4584415038853  44.57942869145138]
[ 4.628901491283342  6.155009451795843 13.432724217601345 
9.365141776937618  44.57942869145138 11.108851081705541]
[  19.7282083595883  3.588909892879645 12.781684520058807  8.381474480151232  38.70544003360641  4.628901491283342]
[ 3.588909892879647  335.6612476370509 27.203654694391926 26.892627599243863  95.61386263390044  6.155009451795841]
[12.781684520058807 27.203654694391926  379.6266330602812  144.8374291115312 384.48762864944314 13.432724217601342]
[ 8.381474480151232 26.892627599243863  144.8374291115312 101.00042848141132 207.84909892879648  9.365141776937618]
[ 38.70544003360641  95.61386263390044 384.48762864944314 207.84909892879648 1536.4584415038853  44.57942869145138]
[ 4.628901491283342  6.155009451795843 13.432724217601345  9.365141776937618  44.57942869145138 11.108851081705541]
T = S.eigenvectors_right() 
       
for t in T: print 'eigenvector of length',t[1][0].norm() print 'principal component of size',t[0] 
       
eigenvector of length 1.0
principal component of size 1699.2947474
eigenvector of length 1.0
principal component of size 8.48911950135
eigenvector of length 1.0
principal component of size 19.5364439354
eigenvector of length 1.0
principal component of size 37.6544565
eigenvector of length 1.0
principal component of size 290.077921102
eigenvector of length 1.0
principal component of size 328.531121681
eigenvector of length 1.0
principal component of size 1699.2947474
eigenvector of length 1.0
principal component of size 8.48911950135
eigenvector of length 1.0
principal component of size 19.5364439354
eigenvector of length 1.0
principal component of size 37.6544565
eigenvector of length 1.0
principal component of size 290.077921102
eigenvector of length 1.0
principal component of size 328.531121681

This just plots the first three coordinates of each point.

points(B.columns()) 
       
Eigs = sorted(T,key=lambda t:t[0],reverse=True) 
       
Eigs 
       
[(1699.294747404675,
  [(0.024891506624188802, 0.07494852053934697, 0.29262379779350217,
0.15039385641720174, 0.9405917591209678, 0.028342119848981332)],
  1),
 (328.5311216810543,
  [(0.0019856117881852328, 0.9955718259783792, -0.002312648556924705,
0.03885283791587622, -0.08512872452228601, 0.008426480145710282)],
  1),
 (290.0779211015012,
  [(0.004944422342956438, -0.03818102812689184, 0.8894781397168148,
0.31908353854771215, -0.32485587780841585, 0.000868807214962941)],
  1),
 (37.65445649996156,
  [(0.10853354340327426, -0.041407655745685126, -0.34927418394187787,
0.9245766540682435, -0.04143206450132864, 0.0891897708083215)],
  1),
 (19.53644393537915,
  [(-0.9437659587213554, 0.0013293000531988486, -0.032108670183043565,
0.12888110266843, 0.02334724876722217, -0.301857000530346)],
  1),
 (8.489119501351412,
  [(0.3112528040782755, 0.006731024084629334, -0.013083627871727066,
0.05104422095259068, 0.015722327894270385, -0.9487109675042368)],
  1)]
[(1699.294747404675,
  [(0.024891506624188802, 0.07494852053934697, 0.29262379779350217, 0.15039385641720174, 0.9405917591209678, 0.028342119848981332)],
  1),
 (328.5311216810543,
  [(0.0019856117881852328, 0.9955718259783792, -0.002312648556924705, 0.03885283791587622, -0.08512872452228601, 0.008426480145710282)],
  1),
 (290.0779211015012,
  [(0.004944422342956438, -0.03818102812689184, 0.8894781397168148, 0.31908353854771215, -0.32485587780841585, 0.000868807214962941)],
  1),
 (37.65445649996156,
  [(0.10853354340327426, -0.041407655745685126, -0.34927418394187787, 0.9245766540682435, -0.04143206450132864, 0.0891897708083215)],
  1),
 (19.53644393537915,
  [(-0.9437659587213554, 0.0013293000531988486, -0.032108670183043565, 0.12888110266843, 0.02334724876722217, -0.301857000530346)],
  1),
 (8.489119501351412,
  [(0.3112528040782755, 0.006731024084629334, -0.013083627871727066, 0.05104422095259068, 0.015722327894270385, -0.9487109675042368)],
  1)]
Basis = [t[1][0] for t in Eigs] 
       
Basis 
       
[(0.024891506624188802, 0.07494852053934697, 0.29262379779350217,
0.15039385641720174, 0.9405917591209678, 0.028342119848981332),
 (0.0019856117881852328, 0.9955718259783792, -0.002312648556924705,
0.03885283791587622, -0.08512872452228601, 0.008426480145710282),
 (0.004944422342956438, -0.03818102812689184, 0.8894781397168148,
0.31908353854771215, -0.32485587780841585, 0.000868807214962941),
 (0.10853354340327426, -0.041407655745685126, -0.34927418394187787,
0.9245766540682435, -0.04143206450132864, 0.0891897708083215),
 (-0.9437659587213554, 0.0013293000531988486, -0.032108670183043565,
0.12888110266843, 0.02334724876722217, -0.301857000530346),
 (0.3112528040782755, 0.006731024084629334, -0.013083627871727066,
0.05104422095259068, 0.015722327894270385, -0.9487109675042368)]
[(0.024891506624188802, 0.07494852053934697, 0.29262379779350217, 0.15039385641720174, 0.9405917591209678, 0.028342119848981332),
 (0.0019856117881852328, 0.9955718259783792, -0.002312648556924705, 0.03885283791587622, -0.08512872452228601, 0.008426480145710282),
 (0.004944422342956438, -0.03818102812689184, 0.8894781397168148, 0.31908353854771215, -0.32485587780841585, 0.000868807214962941),
 (0.10853354340327426, -0.041407655745685126, -0.34927418394187787, 0.9245766540682435, -0.04143206450132864, 0.0891897708083215),
 (-0.9437659587213554, 0.0013293000531988486, -0.032108670183043565, 0.12888110266843, 0.02334724876722217, -0.301857000530346),
 (0.3112528040782755, 0.006731024084629334, -0.013083627871727066, 0.05104422095259068, 0.015722327894270385, -0.9487109675042368)]
P = matrix(Basis).transpose() 
       
       
[ 0.024891506624188802 0.0019856117881852328  0.004944422342956438  
0.10853354340327426   -0.9437659587213554    0.3112528040782755]
[  0.07494852053934697    0.9955718259783792  -0.03818102812689184
-0.041407655745685126 0.0013293000531988486  0.006731024084629334]
[  0.29262379779350217 -0.002312648556924705    0.8894781397168148 
-0.34927418394187787 -0.032108670183043565 -0.013083627871727066]
[  0.15039385641720174   0.03885283791587622   0.31908353854771215   
0.9245766540682435      0.12888110266843   0.05104422095259068]
[   0.9405917591209678  -0.08512872452228601  -0.32485587780841585 
-0.04143206450132864   0.02334724876722217  0.015722327894270385]
[ 0.028342119848981332  0.008426480145710282  0.000868807214962941   
0.0891897708083215    -0.301857000530346   -0.9487109675042368]
[ 0.024891506624188802 0.0019856117881852328  0.004944422342956438   0.10853354340327426   -0.9437659587213554    0.3112528040782755]
[  0.07494852053934697    0.9955718259783792  -0.03818102812689184 -0.041407655745685126 0.0013293000531988486  0.006731024084629334]
[  0.29262379779350217 -0.002312648556924705    0.8894781397168148  -0.34927418394187787 -0.032108670183043565 -0.013083627871727066]
[  0.15039385641720174   0.03885283791587622   0.31908353854771215    0.9245766540682435      0.12888110266843   0.05104422095259068]
[   0.9405917591209678  -0.08512872452228601  -0.32485587780841585  -0.04143206450132864   0.02334724876722217  0.015722327894270385]
[ 0.028342119848981332  0.008426480145710282  0.000868807214962941    0.0891897708083215    -0.301857000530346   -0.9487109675042368]
Y = P.transpose()*B 
       
points([y[:2] for y in Y.columns()]) 
       

So far, so good... but wait a minute!  Recall "What do we actually know about this data set? Right: absolutely nothing!"  So let's use a tool more suited to this.  Here we already do one useful thing - scale the data; but see also here - data analysis is a lot more than just linear algebra or programming alone!

%r liver <- read.csv("/home/notebook/sage_notebook.sagenb/home/kcrisman/130/data/liver.csv") pc <- prcomp( liver, scale=TRUE) summary(pc) 
       
Importance of components:
                          PC1    PC2    PC3    PC4     PC5     PC6    
PC7
Standard deviation     1.5847 1.0946 0.9986 0.9476 0.82032 0.70636
0.47270
Proportion of Variance 0.3588 0.1712 0.1425 0.1283 0.09613 0.07128
0.03192
Cumulative Proportion  0.3588 0.5299 0.6724 0.8007 0.89680 0.96808
1.00000
Importance of components:
                          PC1    PC2    PC3    PC4     PC5     PC6     PC7
Standard deviation     1.5847 1.0946 0.9986 0.9476 0.82032 0.70636 0.47270
Proportion of Variance 0.3588 0.1712 0.1425 0.1283 0.09613 0.07128 0.03192
Cumulative Proportion  0.3588 0.5299 0.6724 0.8007 0.89680 0.96808 1.00000

R graphics don't work so well in Gordon's Sage server, but here it is: