Symbolic Matrix Computation

This example shows how to perform simple matrix computations using Symbolic Math Toolbox™.

Generate a possibly familiar test matrix, the 5-by-5 Hilbert matrix.

H = sym(hilb(5))

H = 
 $(\begin{array}{c} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} \\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} \\ \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} \end{array})$

The determinant is very small.

d = det(H)

d = 
 $\frac{1}{266716800000}$

The elements of the inverse are integers.

X = inv(H)

X = 
 $(\begin{array}{c} 25 & - 300 & 1050 & - 1400 & 630 \\ - 300 & 4800 & - 18900 & 26880 & - 12600 \\ 1050 & - 18900 & 79380 & - 117600 & 56700 \\ - 1400 & 26880 & - 117600 & 179200 & - 88200 \\ 630 & - 12600 & 56700 & - 88200 & 44100 \end{array})$

Verify that the inverse is correct.

I = X*H

I = 
 $(\begin{array}{c} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{array})$

Find the characteristic polynomial.

syms x; p = charpoly(H,x)

p = 
 $x^{5} - \frac{563 x^{4}}{315} + \frac{735781 x^{3}}{2116800} - \frac{852401 x^{2}}{222264000} + \frac{61501 x}{53343360000} - \frac{1}{266716800000}$

Try to factor the characteristic polynomial.

factor(p)

ans = 
 $(\begin{array}{c} \frac{1}{266716800000} & 266716800000 x^{5} - 476703360000 x^{4} + 92708406000 x^{3} - 1022881200 x^{2} + 307505 x - 1 \end{array})$

The result indicates that the characteristic polynomial cannot be factored over the rational numbers.

Compute the 50 digit numerical approximations to the eigenvalues.

digits(50) 
e = eig(vpa(H))

e = 
 $(\begin{array}{c} 1.5670506910982307955330110055207246339493152522334 \\ 0.20853421861101333590500251006882005503858202260343 \\ 0.01140749162341980655945145886658934504234843052664 \\ 0.00030589804015119172687949784069272282565614514909247 \\ 0.0000032879287721718629571150047605447313997367890230746 \end{array})$

Create a generalized Hilbert matrix involving a free variable, $t$ .

t = sym('t'); 
[I,J] = meshgrid(1:5); 
H = 1./(I+J-t)

H = 
 $\begin{array}{l} (\begin{array}{c} - \frac{1}{t - 2} & - \frac{1}{t - 3} & σ_{5} & σ_{3} & σ_{1} \\ - \frac{1}{t - 3} & σ_{5} & σ_{3} & σ_{1} & σ_{2} \\ σ_{5} & σ_{3} & σ_{1} & σ_{2} & σ_{4} \\ σ_{3} & σ_{1} & σ_{2} & σ_{4} & - \frac{1}{t - 9} \\ σ_{1} & σ_{2} & σ_{4} & - \frac{1}{t - 9} & - \frac{1}{t - 10} \end{array}) \\ where \\ σ_{1} = - \frac{1}{t - 6} \\ σ_{2} = - \frac{1}{t - 7} \\ σ_{3} = - \frac{1}{t - 5} \\ σ_{4} = - \frac{1}{t - 8} \\ σ_{5} = - \frac{1}{t - 4} \end{array}$

Substituting $t = 1$ retrieves the original Hilbert matrix.

subs(H,t,1)

ans = 
 $(\begin{array}{c} 1 & \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} \\ \frac{1}{2} & \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} \\ \frac{1}{3} & \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} \\ \frac{1}{4} & \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} \\ \frac{1}{5} & \frac{1}{6} & \frac{1}{7} & \frac{1}{8} & \frac{1}{9} \end{array})$

The reciprocal of the determinant is a polynomial in $t$ .

d = 1/det(H)

d = 
 $- \frac{(t - 2) {(t - 3)}^{2} {(t - 4)}^{3} {(t - 5)}^{4} {(t - 6)}^{5} {(t - 7)}^{4} {(t - 8)}^{3} {(t - 9)}^{2} (t - 10)}{82944}$

d = expand(d)

d = 
 $- \frac{t^{25}}{82944} + \frac{25 t^{24}}{13824} - \frac{5375 t^{23}}{41472} + \frac{40825 t^{22}}{6912} - \frac{15940015 t^{21}}{82944} + \frac{21896665 t^{20}}{4608} - \frac{240519875 t^{19}}{2592} + \frac{1268467075 t^{18}}{864} - \frac{1588946776255 t^{17}}{82944} + \frac{2885896606895 t^{16}}{13824} - \frac{79493630114675 t^{15}}{41472} + \frac{34372691161375 t^{14}}{2304} - \frac{8194259295156385 t^{13}}{82944} + \frac{7707965729450845 t^{12}}{13824} - \frac{55608098247105175 t^{11}}{20736} + \frac{37909434298793825 t^{10}}{3456} - \frac{197019820623693025 t^{9}}{5184} + \frac{10640296363350955 t^{8}}{96} - \frac{38821472549340925 t^{7}}{144} + \frac{12958201048605475 t^{6}}{24} - \frac{1748754621252377 t^{5}}{2} + 1115685328012530 t^{4} - 1078920141906600 t^{3} + 742618453752000 t^{2} - 323874210240000 t + 67212633600000$

The elements of the inverse are also polynomials in $t$ .

X = inv(H)

X = 
 $\begin{array}{l} (\begin{array}{c} - \frac{(t - 2) (t - 3) (t - 4) (t - 5) (t - 6) σ_{4}}{576} & \frac{(t - 3) (t - 4) (t - 5) (t - 6) (t - 7) (t^{4} - 17 t^{3} + 104 t^{2} - 268 t + 240)}{144} & - \frac{(t - 4) (t - 5) (t - 6) (t - 7) (t - 8) (t^{4} - 16 t^{3} + 91 t^{2} - 216 t + 180)}{96} & \frac{(t - 5) (t - 6) (t - 7) (t - 8) (t - 9) (t^{4} - 15 t^{3} + 80 t^{2} - 180 t + 144)}{144} & - \frac{(t - 6) (t - 7) (t - 8) (t - 9) (t - 10) (t^{4} - 14 t^{3} + 71 t^{2} - 154 t + 120)}{576} \\ \frac{(t - 2) (t - 3) (t - 4) (t - 5) (t - 6) σ_{3}}{144} & - \frac{(t - 3) (t - 4) (t - 5) (t - 6) (t - 7) (t^{4} - 21 t^{3} + 161 t^{2} - 531 t + 630)}{36} & \frac{(t - 4) (t - 5) (t - 6) (t - 7) (t - 8) (t^{4} - 20 t^{3} + 145 t^{2} - 450 t + 504)}{24} & - \frac{(t - 5) (t - 6) (t - 7) (t - 8) (t - 9) (t^{4} - 19 t^{3} + 131 t^{2} - 389 t + 420)}{36} & \frac{(t - 6) (t - 7) (t - 8) (t - 9) (t - 10) σ_{4}}{144} \\ - \frac{(t - 2) (t - 3) (t - 4) (t - 5) (t - 6) σ_{2}}{96} & \frac{(t - 3) (t - 4) (t - 5) (t - 6) (t - 7) (t^{4} - 25 t^{3} + 230 t^{2} - 920 t + 1344)}{24} & - \frac{(t - 4) (t - 5) (t - 6) (t - 7) (t - 8) (t^{4} - 24 t^{3} + 211 t^{2} - 804 t + 1120)}{16} & \frac{(t - 5) (t - 6) (t - 7) (t - 8) (t - 9) (t^{4} - 23 t^{3} + 194 t^{2} - 712 t + 960)}{24} & - \frac{(t - 6) (t - 7) (t - 8) (t - 9) (t - 10) σ_{3}}{96} \\ \frac{(t - 2) (t - 3) (t - 4) (t - 5) (t - 6) σ_{1}}{144} & - \frac{(t - 3) (t - 4) (t - 5) (t - 6) (t - 7) (t^{4} - 29 t^{3} + 311 t^{2} - 1459 t + 2520)}{36} & \frac{(t - 4) (t - 5) (t - 6) (t - 7) (t - 8) (t^{4} - 28 t^{3} + 289 t^{2} - 1302 t + 2160)}{24} & - \frac{(t - 5) (t - 6) (t - 7) (t - 8) (t - 9) (t^{4} - 27 t^{3} + 269 t^{2} - 1173 t + 1890)}{36} & \frac{(t - 6) (t - 7) (t - 8) (t - 9) (t - 10) σ_{2}}{144} \\ - \frac{(t - 2) (t - 3) (t - 4) (t - 5) (t - 6) (t^{4} - 34 t^{3} + 431 t^{2} - 2414 t + 5040)}{576} & \frac{(t - 3) (t - 4) (t - 5) (t - 6) (t - 7) (t^{4} - 33 t^{3} + 404 t^{2} - 2172 t + 4320)}{144} & - \frac{(t - 4) (t - 5) (t - 6) (t - 7) (t - 8) (t^{4} - 32 t^{3} + 379 t^{2} - 1968 t + 3780)}{96} & \frac{(t - 5) (t - 6) (t - 7) (t - 8) (t - 9) (t^{4} - 31 t^{3} + 356 t^{2} - 1796 t + 3360)}{144} & - \frac{(t - 6) (t - 7) (t - 8) (t - 9) (t - 10) σ_{1}}{576} \end{array}) \\ where \\ σ_{1} = t^{4} - 30 t^{3} + 335 t^{2} - 1650 t + 3024 \\ σ_{2} = t^{4} - 26 t^{3} + 251 t^{2} - 1066 t + 1680 \\ σ_{3} = t^{4} - 22 t^{3} + 179 t^{2} - 638 t + 840 \\ σ_{4} = t^{4} - 18 t^{3} + 119 t^{2} - 342 t + 360 \end{array}$

Substituting $t = 1$ generates the Hilbert inverse.

X = subs(X,t,'1')

X = 
 $(\begin{array}{c} 25 & - 300 & 1050 & - 1400 & 630 \\ - 300 & 4800 & - 18900 & 26880 & - 12600 \\ 1050 & - 18900 & 79380 & - 117600 & 56700 \\ - 1400 & 26880 & - 117600 & 179200 & - 88200 \\ 630 & - 12600 & 56700 & - 88200 & 44100 \end{array})$

X = double(X)

X = 5×5

          25        -300        1050       -1400         630
        -300        4800      -18900       26880      -12600
        1050      -18900       79380     -117600       56700
       -1400       26880     -117600      179200      -88200
         630      -12600       56700      -88200       44100

Investigate a different example.

A = sym(gallery(5))

A = 
 $(\begin{array}{c} - 9 & 11 & - 21 & 63 & - 252 \\ 70 & - 69 & 141 & - 421 & 1684 \\ - 575 & 575 & - 1149 & 3451 & - 13801 \\ 3891 & - 3891 & 7782 & - 23345 & 93365 \\ 1024 & - 1024 & 2048 & - 6144 & 24572 \end{array})$

This matrix is "nilpotent". Its fifth power is the zero matrix.

A^5

ans = 
 $(\begin{array}{c} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \end{array})$

Because this matrix is nilpotent, its characteristic polynomial is very simple.

p = charpoly(A,'lambda')

p =  $λ^{5}$

You should now be able to compute the matrix eigenvalues in your head. They are the zeros of the equation lambda^5 = 0.

Symbolic computation can find the eigenvalues exactly.

lambda = eig(A)

lambda = 
 $(\begin{array}{c} 0 \\ 0 \\ 0 \\ 0 \\ 0 \end{array})$

Numeric computation involves roundoff error and finds the zeros of an equation that is something like lambda^5 = eps*norm(A) So the computed eigenvalues are roughly lambda = (eps*norm(A))^(1/5) Here are the eigenvalues, computed by the Symbolic Toolbox using 16 digit floating point arithmetic. It is not obvious that they should all be zero.

digits(16) 
lambda = eig(vpa(A))

lambda = 
 $(\begin{array}{c} 0.0005617448486395847 \\ 0.0001737348850386136 - 0.0005342985684139864 i \\ 0.0001737348850386136 + 0.0005342985684139864 i \\ - 0.000454607309358406 + 0.0003303865815979566 i \\ - 0.000454607309358406 - 0.0003303865815979566 i \end{array})$

This matrix is also "defective". It is not similar to a diagonal matrix. Its Jordan Canonical Form is not diagonal.

J = jordan(A)

J = 
 $(\begin{array}{c} 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 \end{array})$

The matrix exponential, expm(t*A), is usually expressed in terms of scalar exponentials involving the eigenvalues, exp(lambda(i)*t). But for this matrix, the elements of expm(t*A) are all polynomials in t.

t = sym('t'); 
E = simplify(expm(t*A))

E = 
 $(\begin{array}{c} - \frac{2 t^{3}}{3} + \frac{11 t^{2}}{2} - 9 t + 1 & \frac{t (4 t^{2} - 27 t + 33)}{3} & - \frac{t (20 t^{2} - 117 t + 126)}{6} & \frac{t (32 t^{2} - 174 t + 189)}{3} & - \frac{t (85 t^{2} - 464 t + 504)}{2} \\ - \frac{t (7 t^{3} - 81 t^{2} + 230 t - 140)}{2} & 7 t^{4} - 67 t^{3} + \frac{301 t^{2}}{2} - 69 t + 1 & - \frac{t (35 t^{3} - 293 t^{2} + 598 t - 282)}{2} & \frac{t (112 t^{3} - 876 t^{2} + 1799 t - 842)}{2} & - \frac{t (1785 t^{3} - 14012 t^{2} + 28776 t - 13472)}{8} \\ \frac{t (142 t^{3} - 1710 t^{2} + 5151 t - 3450)}{6} & - \frac{t (142 t^{3} - 1426 t^{2} + 3438 t - 1725)}{3} & \frac{355 t^{4}}{3} - \frac{3139 t^{3}}{3} + \frac{4585 t^{2}}{2} - 1149 t + 1 & - \frac{t (1136 t^{3} - 9420 t^{2} + 20625 t - 10353)}{3} & \frac{t (18105 t^{3} - 150646 t^{2} + 329952 t - 165612)}{12} \\ - \frac{t (973 t^{3} - 11675 t^{2} + 35022 t - 23346)}{6} & \frac{t (1946 t^{3} - 19458 t^{2} + 46695 t - 23346)}{6} & - \frac{t (4865 t^{3} - 42807 t^{2} + 93390 t - 46692)}{6} & \frac{7784 t^{4}}{3} - \frac{64210 t^{3}}{3} + \frac{93391 t^{2}}{2} - 23345 t + 1 & - \frac{t (248115 t^{3} - 2053748 t^{2} + 4482036 t - 2240760)}{24} \\ - \frac{128 t (t^{3} - 12 t^{2} + 36 t - 24)}{3} & \frac{256 t (t^{3} - 10 t^{2} + 24 t - 12)}{3} & - \frac{128 t (5 t^{3} - 44 t^{2} + 96 t - 48)}{3} & \frac{512 t (4 t^{3} - 33 t^{2} + 72 t - 36)}{3} & - 2720 t^{4} + \frac{67552 t^{3}}{3} - 49144 t^{2} + 24572 t + 1 \end{array})$

By the way, the function "exp" computes element-by-element exponentials.

X = exp(t*A)

X = 
 $(\begin{array}{c} e^{- 9 t} & e^{11 t} & e^{- 21 t} & e^{63 t} & e^{- 252 t} \\ e^{70 t} & e^{- 69 t} & e^{141 t} & e^{- 421 t} & e^{1684 t} \\ e^{- 575 t} & e^{575 t} & e^{- 1149 t} & e^{3451 t} & e^{- 13801 t} \\ e^{3891 t} & e^{- 3891 t} & e^{7782 t} & e^{- 23345 t} & e^{93365 t} \\ e^{1024 t} & e^{- 1024 t} & e^{2048 t} & e^{- 6144 t} & e^{24572 t} \end{array})$