# MAT 232 Chapter Five

## 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)
 \newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & 1 & 3 \\ 1 & 3 & 1 \\ 3 & 1 & 1 \end{array}\right)
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.

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: