Lay Chapter 7,
Symmetric Matrices and Quadratic Forms

Diagonalize a Matrix

In[1]:=

Clear[a] ; 

a = ( {{3, 1}, {1, 3}} ) ;

Calculate the eigenvalues of a.

In[3]:=

evals = Eigenvalues[a] ;

Λ = DiagonalMatrix[evals] ;

%//MatrixForm

Out[5]//MatrixForm=

( {{4, 0}, {0, 2}} )

Calculate the eigenvectors of a.

In[6]:=

evecs = Eigenvectors[a] ;

v = Transpose[evecs] ;

%//MatrixForm

Out[8]//MatrixForm=

( {{1, -1}, {1, 1}} )

Normalize the columns of v

In[9]:=

v1 = Transpose[v][[1]] ;

v2 = Transpose[v][[2]] ;

p = Transpose[{1/(v1 . v1)^(1/2) v1, 1/(v2 . v2)^(1/2) v2}] ;

%//MatrixForm

Out[12]//MatrixForm=

( {{1/2^(1/2), -1/2^(1/2)}, {1/2^(1/2), 1/2^(1/2)}} )

Check.

In[13]:=

Transpose[p] Inverse[p] ;

a  p . Λ . Transpose[p]

Out[14]=

True

Diagonalize a Matrix

In[15]:=

Clear[a] ; 

a = ({{5, 2, 9, -6}, {2, 5, -6, 9}, {9, -6, 5, 2}, {-6, 9, 2, 5}}) ;

Calculate the eigenvalues of a.

In[17]:=

evals = Eigenvalues[a] ;

Λ = DiagonalMatrix[evals] ;

%//MatrixForm

Out[19]//MatrixForm=

( {{18, 0, 0, 0}, {0, -12, 0, 0}, {0, 0, 10, 0}, {0, 0, 0, 4}} )

Calculate the eigenvectors of a.

In[20]:=

evecs = Eigenvectors[a] ;

v = Transpose[evecs] ;

%//MatrixForm

Out[22]//MatrixForm=

( {{-1, 1, 1, -1}, {1, -1, 1, -1}, {-1, -1, 1, 1}, {1, 1, 1, 1}} )

Normalize the columns of v

In[23]:=

v1 = Transpose[v][[1]] ;

v2 = Transpose[v][[2]] ;

v3 = Transpose[v][[3]] ;

v4 = Transpose[v][[4]] ;

p = Transpose[{1/(v1 . v1)^(1/2) v1, 1/(v2 . v2)^(1/2) v2, 1/(v3 . v3)^(1/2) v3, 1/(v4 . v4)^(1/2) v4}] ;

%//MatrixForm

Out[28]//MatrixForm=

( {{-1/2, 1/2, 1/2, -1/2}, {1/2, -1/2, 1/2, -1/2}, {-1/2, -1/2, 1/2, 1/2}, {1/2, 1/2, 1/2, 1/2}} )

Check.

In[29]:=

Transpose[p] Inverse[p]

p . Λ . Transpose[p] a

Out[29]=

True

Out[30]=

True

Matrix Representation of a Quadratic Form

Symmetric matrices represent quadratic forms in the following way.

In[31]:=

Clear[a, x, x1, x2] 

a = ( {{1, 2}, {2, 3}} ) ;

x = {x1, x2} ; 

x . a . x//Simplify

Out[34]=

x1^2 + 4 x1 x2 + 3 x2^2

Here is the quadratic form generated by a diagonal matrix

In[35]:=

Clear[a, x, x1, x2] 

a = ( {{4, 0}, {0, 3}} ) ;

x = {x1, x2} ; 

x . a . x//Simplify

Out[38]=

4 x1^2 + 3 x2^2

... and another generated by a symmetric matrix which is not a diagonal matrix.

In[39]:=

Clear[a, x, x1, x2] 

a = ( {{3, -2}, {-2, 7}} ) ;

x = {x1, x2} ; 

x . a . x//Simplify

Out[42]=

3 x1^2 - 4 x1 x2 + 7 x2^2

The difference is the "mixed" term.

Matrix Representation of Quadratic Forms

Write the matrix representing a particular quadratic form.

In[43]:=

Clear[a, x, x1, x2, x3] 

a = ({{5, -1/2, 0}, {-1/2, 3, 4}, {0, 4, 2}}) ;

x = {x1, x2, x3} ; 

x . a . x//Simplify

Out[46]=

5 x1^2 - x1 x2 + 3 x2^2 + 8 x2 x3 + 2 x3^2

Value of a Quadratic Form at a Point

Evaluate a particular quadratic form at distinct points.

In[47]:=

Clear[a, x, x1, x2, q] 

a = ( {{1, -4}, {-4, -5}} ) ;

x = {x1, x2} ; 

q[x_] := x . a . x//Simplify

In[51]:=

q[{-3, 1}]

q[{2, -2}]

q[{1, -3}]

Out[51]=

28

Out[52]=

16

Out[53]=

-20

Change of Variable in a Quadratic Form

If x is a variable in ^n, thena change of variable in a quadratic form is an equation of the form.

In[54]:=

(* x = p y *)

where y is a variable in ^n and p is an invertible n x n matrix.

Diagonalize the following quadratic form.

In[55]:=

Clear[a, x, x1, x2, q, evals, evecs, d, p] 

a = ( {{1, -4}, {-4, -5}} ) ;

x = {x1, x2} ; 

qx[x_] := x . a . x//Simplify

First, calculate the corresponding diagonal matrix.

In[59]:=

evals = Eigenvalues[a] ;

d = DiagonalMatrix[evals] ;

%//MatrixForm

Out[61]//MatrixForm=

( {{-7, 0}, {0, 3}} )

Then the columns of the matrix p are the orthonormal eigenvectors of a, in the proper order.

In[62]:=

evecs = Eigenvectors[a] ;

d1 = evecs[[1]] ;

d2 = evecs[[2]] ; 

p = Transpose[{1/(d1 . d1)^(1/2) d1, 1/(d2 . d2)^(1/2) d2}] ;

%//MatrixForm

Out[66]//MatrixForm=

( {{1/5^(1/2), -2/5^(1/2)}, {2/5^(1/2), 1/5^(1/2)}} )

Check the diagonalization.

In[67]:=

ap . d . Inverse[p]

Out[67]=

True

Check that the two quadratic forms agree for x = {2,-2}.
First, here is qx applied to x = {2,-2}.

In[68]:=

qx[{2, -2}]

Out[68]=

16

Now x = p.y implies that y = Inverse[p].x
Here is qy applied to y = Inverse[p].x

In[69]:=

qy[y_] := y . d . y

qy[Inverse[p] . {2, -2}]

Out[70]=

16

Colorful Vectors and Curves
Contains Mathematica Code Only -- No Solved Exercises

Principal Axes

Let's recast a quadratic form generated by a symmetric 2 x 2 matrix as a function of two variables.
We use the matrix a from the previous example.

In[97]:=

Clear[qa, qd, a, d, x1, x2, y1, y2] 

a = ( {{1, -4}, {-4, -5}} ) ; 

qa[x1_, x2_] := {x1, x2} . a . {x1, x2}

In[100]:=

qaSurface = Plot3D[qa[x1, x2], {x1, -3, 3}, {x2, -3, 3}, AxesLabel {"x", "y", "z"}] ;

[Graphics:HTMLFiles/7_quadratic_forms_133.gif]

Now do the same thing for the matrix d.

In[101]:=

d = ( {{-7, 0}, {0, 3}} ) ; 

qd[y1_, y2_] := {y1, y2} . d . {y1, y2}

In[103]:=

qdSurface = Plot3D[qd[y1, y2], {y1, -3, 3}, {y2, -3, 3}, AxesLabel {"x", "y", "z"}] ;

[Graphics:HTMLFiles/7_quadratic_forms_137.gif]

The second surface s obtained from the first by a rotation.
The y1- and y2-axes in the second figure correspond to the principal axes in the first.
The principal axes are the directions indicated by the columns of the matrix p used in describing the change of variables: x = p y.

In[104]:=

origin = {0, 0, 0} ;

p1 = {Transpose[p][[1]], 0}//Flatten

p2 = {Transpose[p][[2]], 0}//Flatten

Out[105]=

{1/5^(1/2), 2/5^(1/2), 0}

Out[106]=

{-2/5^(1/2), 1/5^(1/2), 0}

Let's see those two vectors. We will exaggerate their lengths in order to see them better.
Note the orientation of the standard coordinate axes.

In[107]:=

showColorful3DVectors[{{origin, 5 p1}, {origin, 5 p2}}, Red, Tomato, PlotRangeAll] ;

[Graphics:HTMLFiles/7_quadratic_forms_144.gif]

In this image as well, the lengths of the vectors indicating the principal axes are exagerated in order to see them better.
One of the vectors is mostly underneath the surface.

In[108]:=

principalAxes = showColorful3DVectors[{{origin, 5 p1}, {origin, 5 p2}}, Red, Orchid, PlotRangeAll, DisplayFunctionIdentity] ; 

Show[qaSurface, principalAxes, PlotRangeAll, AxesFalse, BoxedFalse, DisplayFunction$DisplayFunction] ;

[Graphics:HTMLFiles/7_quadratic_forms_147.gif]

Using Contour Plots to Investigate Quadratic Forms and their Principal Axes

Let's see a contour plot of the quadratic form qa.

In[110]:=

Clear[qa, qd, a, d, x1, x2, y1, y2] 

a = ( {{1, -4}, {-4, -5}} ) ; 

qa[x1_, x2_] := {x1, x2} . a . {x1, x2} 

contours = ContourPlot[qa[x1, x2], {x1, -3, 3}, {x2, -3, 3}] ;

[Graphics:HTMLFiles/7_quadratic_forms_152.gif]

... and a contour plot of the quadratic form qd.
It seems clear that the second is obtained by a rotation of the first.

In[114]:=

d = ( {{-7, 0}, {0, 3}} ) ; 

qd[y1_, y2_] := {y1, y2} . d . {y1, y2} 

ContourPlot[qd[y1, y2], {y1, -3, 3}, {y2, -3, 3}] ;

[Graphics:HTMLFiles/7_quadratic_forms_156.gif]

Now let's see those principal axes superimposed on the contour plot of qa.

In[117]:=

origin = {0, 0} ;

pa1 = {origin, Transpose[p][[1]]} ;

pa2 = {origin, Transpose[p][[2]]} ; 

vecs = showColorfulVectors[{pa1, pa2}, ManganeseBlue, DisplayFunctionIdentity] ; 

Show[contours, vecs, DisplayFunction$DisplayFunction] ;

[Graphics:HTMLFiles/7_quadratic_forms_162.gif]

Classification of Quadratic Forms

Positive definite

In[122]:=

Clear[a, qa, qd] 

a = ( {{2, 0}, {0, 3}} ) ; 

qa[x1_, x2_] := {x1, x2} . a . {x1, x2} 

qaSurface = Plot3D[qa[x1, x2], {x1, -3, 3}, {x2, -3, 3}, AxesLabel {"x", "y", "z"}] ;

[Graphics:HTMLFiles/7_quadratic_forms_167.gif]

Negative definite

In[126]:=

Clear[a, qa, qd] 

a = ( {{-2, 0}, {0, -3}} ) ; 

qa[x1_, x2_] := {x1, x2} . a . {x1, x2} 

qaSurface = Plot3D[qa[x1, x2], {x1, -3, 3}, {x2, -3, 3}, AxesLabel {"x", "y", "z"}] ;

[Graphics:HTMLFiles/7_quadratic_forms_172.gif]

Indefinite

In[130]:=

Clear[a, qa, qd] 

a = ( {{2, 0}, {0, -3}} ) ; 

qa[x1_, x2_] := {x1, x2} . a . {x1, x2} 

qaSurface = Plot3D[qa[x1, x2], {x1, -3, 3}, {x2, -3, 3}, AxesLabel {"x", "y", "z"}] ;

[Graphics:HTMLFiles/7_quadratic_forms_177.gif]

Positive Semidefinite

In[134]:=

Clear[a, qa, qd] 

a = ( {{2, 0}, {0, 0}} ) ; 

qa[x1_, x2_] := {x1, x2} . a . {x1, x2} 

qaSurface = Plot3D[qa[x1, x2], {x1, -3, 3}, {x2, -3, 3}, AxesLabel {"x", "y", "z"}] ;

[Graphics:HTMLFiles/7_quadratic_forms_182.gif]

Diagonalization a Quadratic Form

Diagonalize the following quadratic form.

In[138]:=

Clear[a, x, x1, x2, x3, x4, q, evals, evecs, d, p] 

a = ({{1, 9/2, 0, -6}, {9/2, 1, 6, 0}, {0, 6, 1, 9/2}, {-6, 0, 9/2, 1}}) ;

x = {x1, x2, x3, x4} ; 

qx[x_] = x . a . x//Simplify

Out[141]=

x1^2 + 9 x1 x2 + x2^2 + 12 x2 x3 + x3^2 - 12 x1 x4 + 9 x3 x4 + x4^2

First, calculate the corresponding diagonal matrix.

In[142]:=

evals = Eigenvalues[a] ;

d = DiagonalMatrix[evals] ;

%//MatrixForm

Out[144]//MatrixForm=

( {{17/2, 0, 0, 0}, {0, 17/2, 0, 0}, {0, 0, -13/2, 0}, {0, 0, 0, -13/2}} )

Use d to construct the equivalent (diagonal) quadratic form.
This quadratic form is Lay 7.2.17.

In[145]:=

qd[y1_, y2_, y3_, y4_] = {y1, y2, y3, y4} . d . {y1, y2, y3, y4}

Out[145]=

(17 y1^2)/2 + (17 y2^2)/2 - (13 y3^2)/2 - (13 y4^2)/2

Now relate the two quadratic forms by calculating the change of basis matrix p such that x = p y.
The columns of the matrix p are the orthonormal eigenvectors of a, in the proper order.

In[146]:=

evecs = Eigenvectors[a] ;

d1 = evecs[[1]] ;

d2 = evecs[[2]] ;

d3 = evecs[[3]] ;

d4 = evecs[[4]] ; 

p = Transpose[{1/(d1 . d1)^(1/2) d1, 1/(d2 . d2)^(1/2) d2, 1/(d3 . d3)^(1/2) d3, 1/(d4 . d4)^(1/2) d4}] ;

%//MatrixForm

Out[152]//MatrixForm=

( {{-1/2^(1/2), 3/(5 2^(1/2)), 1/2^(1/2), 3/(5 2^(1/2))}, {-3/(5 2^(1/2)), 1/2^(1/2) ... /2^(1/2)}, {0, (2 2^(1/2))/5, 0, (2 2^(1/2))/5}, {(2 2^(1/2))/5, 0, (2 2^(1/2))/5, 0}} )

Check the diagonalization.

In[153]:=

ap . d . Inverse[p]

Out[153]=

True

Constrained Optimization

Consider the following quadratic form:

In[154]:=

Clear[a, x, x1, x2, x3, evals, evecs, d, p] 

a = ({{3, 2, 1}, {2, 3, 1}, {1, 1, 4}}) ;

x = {x1, x2, x3} ; 

x . a . x//Simplify

Out[157]=

3 x1^2 + 3 x2^2 + 2 x2 x3 + 4 x3^2 + 2 x1 (2 x2 + x3)

Calculate the corresponding diagonal matrix.

In[158]:=

evals = Eigenvalues[a] ;

d = DiagonalMatrix[evals] ;

%//MatrixForm

Out[160]//MatrixForm=

( {{6, 0, 0}, {0, 3, 0}, {0, 0, 1}} )

The columns of the matrix p are the orthonormal eigenvectors of a, taken in the proper order.

In[161]:=

evecs = Eigenvectors[a] ;

d1 = evecs[[1]] ;

d2 = evecs[[2]] ;

d3 = evecs[[3]] ; 

p = Transpose[{1/(d1 . d1)^(1/2) d1, 1/(d2 . d2)^(1/2) d2, 1/(d3 . d3)^(1/2) d3}] ;

%//MatrixForm

Out[166]//MatrixForm=

( {{1/3^(1/2), -1/6^(1/2), -1/2^(1/2)}, {1/3^(1/2), -1/6^(1/2), 1/2^(1/2)}, {1/3^(1/2), 2/3^(1/2), 0}} )

Check the diagonalization.

In[167]:=

ap . d . Inverse[p]

Out[167]=

True

Constrained Optimization

Consider the following quadratic form.

In[168]:=

Clear[a, x, x1, x2, evals, evecs, d, p] 

a = ( {{3, 1}, {1, 3}} ) ;

x = {x1, x2} ; 

x . a . x//Simplify

Out[171]=

3 x1^2 + 2 x1 x2 + 3 x2^2

Calculate the corresponding diagonal matrix.

In[172]:=

evals = Eigenvalues[a] ;

d = DiagonalMatrix[evals] ;

%//MatrixForm

Out[174]//MatrixForm=

( {{4, 0}, {0, 2}} )

Then the columns of the matrix p are the orthonormal eigenvectors of a, taken in the proper order.

In[175]:=

evecs = Eigenvectors[a] ;

d1 = evecs[[1]] ;

d2 = evecs[[2]] ; 

p = Transpose[{1/(d1 . d1)^(1/2) d1, 1/(d2 . d2)^(1/2) d2}] ;

%//MatrixForm

Out[179]//MatrixForm=

( {{1/2^(1/2), -1/2^(1/2)}, {1/2^(1/2), 1/2^(1/2)}} )

Check the diagonalization.

In[180]:=

ap . d . Inverse[p]

Out[180]=

True

We can graph this surface.
Stretch out the graph for a better view.

In[181]:=

Clear[x1, x2] 

q[x1_, x2_] = If [x1^2 + x2^2≤ 1, {x1, x2} . a . {x1, x2}, -1] ; 

Plot3D[q[x1, x2], {x1, -1.1, 1.1}, {x2, -1.1, 1.1}, PlotRange {0, 5}, MeshFalse, PlotPoints100, AxesLabel {x1, x2, None}] ;

[Graphics:HTMLFiles/7_quadratic_forms_242.gif]

Created by Mathematica  (December 11, 2004) Valid XHTML 1.1!