MGM.10 The Spectral Theorem for Symmetric Matrices and the Holy Grail of Matrix Theory
TUTORIALS

Needs["LinearAlgebra`MatrixManipulation`"];
Needs["LinearAlgebra`Orthogonalization`"];
Needs["Graphics`Colors`"];
Off[General::"spell"];
Off[General::"spell1"];
Off[ParametricPlot3D::"ppcom"];
Off[ParametricPlot::"ppcom"];
Off[Plot::"plnr"];
SetOptions[Limit, Analytic -> True];
$Pre = Chop;

CMView = {2.7, 1.6, 1.2};

Vector Drawer

\(\(\( (*\ \(\(:\)\(Date : \ \ Copyright\ 2001\)\), \ 
      Math\ Everywhere, \ \(\(Inc\)\(.\)\)\ 
*) \)\(
\)\( (*
Original\ Conception\ to\ create\ a\ vector\ drawer\ is\ due\ to\ Bill\ Davis . 
Current\ code\ is\ by\ Bruce\ \(\(Carpenter\)\(.\)\)
*) \)\(
\)\( (*\ \ \(\(\(:\)\(Drawbacks : \ \ Since\ the\ vectors\ are\ drawn\ out\ of\)\)\ 

      \ \ \ \ \ context\), \ 
      the\ shapes\ of\ the\ arrowheads\ are\ only\ 
\ \ \ \ \ correct\ with\ AspectRatio\ set\ to\ \(\(Automatic\)\(.\)\)
*) \)\)\)

BeginPackage["Vector3D`","Graphics`Colors`"];

Axes3D::usage = "Axes3D[a, b] creates a Graphics3D object of cartesian axes with x, y, and z running from -a/3 to a, and with axes labels b units beyond the tips of the axes.  Axes3D[a] is Axes3D[a, a/8].";

Perpend::usage = "Perpend[a] returns a unit vector which is perpendicular to the vector a.";

Arrow::usage = "Arrow[a] produces a 2 or 3 dimensional vector from the origin to a. Arrow[a, Tail -> tail] gives a vector from tail to a.";
	
Vector::usage = "Vector[a, b] gives a 2D or 3D vector from a to b.";

ArrowHead::usage = "ArrowHead[a, vec] produces an arrowhead with its tip placed at the point a, pointing in the direction of the vector vec.";

Tail::usage = "Tail -> point puts the tail of the vector at point.";

Norm::usage = "Norm[vec] returns the Euclidean length of vec.";

Normalize::usage = "Normalize[vec] returns a unit vector in the direction of vec, provided vec is nonzero.";

Aperture::usage = "Aperture is the ratio of the radius of the base to the length of the head of a vector.";

TipSize::usage = "TipSize specifies an absolute size for the inner tip length of the head of a vector.";

TipRatio::usage = "TipRatio is the ratio of the inner length of the arrowhead to the length of the shaft of a vector.";

HeadRatio::usage = "HeadRatio is the ratio of the outer length of the arrowhead to the length of the shaft of a vector.";

HeadSize::usage = "HeadSize specifies an absolute size for the head of a vector.";

ScaleFactor::usage = "ScaleFactor specifies the amount to scale a vector in length and may be either a number or a pure function.";

ZeroVectorPointSize::usage = "ZeroVectorPointSize is the size of the point used to represent a zero vector.";

Options[Vector] = Options[Arrow] = Options[ArrowHead] =
{HeadRatio -> 0.18, TipRatio -> 0.14, Aperture -> 0.3, ScaleFactor -> 1, HeadType -> Polygon, ShaftQ -> True, EdgesQ -> True, VectorColor -> RGBColor[0,0,1], ShaftWidth -> 0.005, ZeroVectorPointSize -> 0.01};
		
Begin["`Private`"];

SetOptions[ParametricPlot,AspectRatio -> Automatic];
SetOptions[Plot,AspectRatio -> Automatic];
SetOptions[Graphics,AspectRatio -> Automatic];

Axes3D[u_,v_] := 
      Graphics3D[{
            {Blue,Line[{{-u/3,0,0},{u,0,0}}]},
                    Text["x",{u + v,0,0}],
            {Blue,Line[{{0,-u/3,0},{0,u,0}}]},
                    Text["y",{0, u + v,0}],
            {Blue,Line[{{0,0,-u/3},{0,0,u}}]},
                    Text["z",{0,0,u + v}]}]; 
Axes3D[u_] := Axes3D[u,u/8];
				
Norm[a_?VectorQ] := N[Sqrt[a.a]];

Normalize[a_?VectorQ] := If[a.a == 0, a, a/Norm[a]];

Perpend[a_?VectorQ] := Normalize[{a[[2]],-a[[1]]}]/;Length[a]==2;

Perpend[a_?VectorQ]:=
    Normalize[If[a[[1]]==0,{1,0,0},{a[[2]],-a[[1]],0}]]/;
      Length[a]==3;
      
Base3D = Table[{Cos[2 k Pi/8.], Sin[2 k Pi/8.]},{k,0,8}]; 

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, sw, ht, hs, hr, tr, ap, sf, vc, zvps, diff = to - from, tip},
   {sq, sw, ht, hs, hr, ts, tr, ap, sf, vc, zvps} = {ShaftQ, ShaftWidth, HeadType, HeadSize, HeadRatio, TipSize, TipRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]; tr = 0.75 hr];
   If[NumberQ[ts], tr = ts/Norm[diff]];
   Graphics[{If[sq,{vc,Thickness[sw],Line[{from,tip - tr diff}]},{}],
      With[{trans = ap hr Norm[diff] Perpend[diff]},
      {vc, Thickness[sw], ht[{tip,
                   (tip - hr diff + trans),
                   (tip - tr diff),
                   (tip - hr diff - trans),tip}]}]}]
] /; Length[from]==Length[to]==2

Vector[from_?VectorQ, to_?VectorQ, opts___] := Module[
   {sq, eq, sw, hs, hr, tr,ap, sf, vc, zvps, diff = to - from, tip},
   {sq, eq, sw, hs, hr, ap, sf, vc, zvps} = {ShaftQ, EdgesQ, ShaftWidth, HeadSize, HeadRatio, Aperture, ScaleFactor, VectorColor, ZeroVectorPointSize} /. {opts} /. Options[Vector];
   If[from==to, Return[Graphics3D[{PointSize[zvps],vc,Point[from]}]]];
   diff = If[NumberQ[N[sf]], sf diff, sf[diff]];
   tip = from + diff;
   If[NumberQ[hs], hr = hs/Norm[diff]];
   tr = hr/2;
   Graphics3D[{If[sq,{vc,Thickness[sw],Line[{from,tip-tr diff}]},{}], If[eq, {}, EdgeForm[]], {SurfaceColor[vc],
      With[{temp = Perpend[diff]},
      Polygon[Append[#,tip]]& /@ Partition[(tip - hr diff + #)&/@((ap hr Norm[diff] Base3D) . {temp,Normalize[Cross[temp,diff]]}),2,1]]}}]
] /; Length[from]==Length[to]==3


Arrow[a_, Tail->b_, opts___] := Vector[b,a+b,opts]

Arrow[a_, opts___] := Vector[Table[0,{Length[a]}],a,opts]

ArrowHead[a_, b_, opts___] := Vector[a-b,a,ShaftQ->False,opts];

End[];

EndPackage[];
Null

Graphics`Colors`GosiaGreen = 
        RGBColor[0,  0.392187,  0];

ThreeAxes[u,v]

ThreeAxes[u_, v_] := 
        Graphics3D[{{Blue, Line[{{-u, 0, 0}, {u, 0, 0}}]}, 
            Text["x", {u + v, 0, 0}], {Blue, 
              Line[{{0, -u, 0}, {0, u, 0}}]}, 
            Text["y", {0, u + v, 0}], {Blue, 
              Line[{{0, 0, 0}, {0, 0, u}}]}, 
            Text["z", {0, 0, u + v}]}];
ThreeAxes[u_] := ThreeAxes[u, (u/8)];
ThreeAxes::"usage" = "ThreeAxes[a,b] makes a standard cartesian \
axis graphics object with x, y, and z running from -a to a, and with axis \
labels b units beyond the tips of the axes.  ThreeAxes[a] is \
ThreeAxes[a,a/8].";

Needs["Statistics`ContinuousDistributions`"];

Clear[matrixplucker];
matrixplucker[expr_, 
        vars_] := Simplify[
          Outer[D, {expr}, 
            vars]][[1]]

T.1) If A is of full rank, then PseudoInverse[A] =  ((((A^t) . A))^(-1)) (A^t)
If A is not of full rank, then the formula fails.

T.1.a.i)  If A is of full rank, then PseudoInverse[A] =  ((((A^t) . A))^(-1)) (A^t)

Here is a random full rank matrix A:

hitdim  =  Random[Integer, {3, 7}];
hangdim  = 
        hitdim  +  Random[Integer, {2, 6}];
A  =  
        Table[Random[Real, {-3, 3}], {i, 1, hangdim}, {j, 1, 
            hitdim}];
MatrixForm[A]

Check:

rank = Length[SingularValues[A][[2]]];
rank  ==  hitdim

Now look at Mathematica's calculation of PseudoInverse[A]:

MatrixForm[PseudoInverse[A]]

And Mathematica's calculation of  ((((A^t) . A))^(-1)) (A^t).

MatrixForm[Inverse[Transpose[A] . A] . Transpose[A]]

They are the same. This is superb evidence that 
               PseudoInverse[A] =  ((((A^t) . A))^(-1)) (A^t).
Explain why this is guaranteed for any and all matrices A of full rank.

Go with any full rank matrix A.
Go with any Y in hangdimD. 
         Xclosest = PseudoInverse[A].Y
means A.Xclosest is closer to Y than any other hit with A. 

And because perpendicular distance is the shortest distance, this means 
Y - A.Xclosest is perpendicular to the subspace of hangdimD consisting of all hits with A. 

In other words  (Y - A.Xclosest).(A.X) = 0 for all X's in hitdimD.

In view of the transpose manipulation, this is the same as 
           ((A^t) .(Y - A.Xclosest)).X = 0 for all X in hitdimD.
           
The only vector in hitdimD perpendicular to all the X's in hitdimD is {0,0,,...,0}.

So      {0,0,. . .,0} =  (A^t) .(Y - A.Xclosest) = (A^t).Y -  (A^t).A.Xclosest)).

This is the same as
            (A^t).Y =  (A^t).A.Xclosest.
Because A is of full rank, (A^t).A is invertible.
 
Hit both sides with   ((((A^t) . A))^(-1)) to get
  ((((A^t) . A))^(-1)).(A^t).Y =  ((((A^t) . A))^(-1)).(A^t).A.Xclosest  = Xclosest. = PseudoInverse[A].Y.
  
  Read across to see that
         ((((A^t) . A))^(-1)).(A^t).Y  = PseudoInverse[A].Y.
Because this happens for any Y in hangdimD, this tells you that
          ((((A^t) . A))^(-1)).(A^t) = PseudoInverse[A].
And you're out of here.

T.1.a.ii)  If A is of not of full rank, then the formula fails

Here is a full rank matrix A:

hitdim  =  Random[Integer, {4, 7}];
hangdim  = 
        hitdim  -  Random[Integer, {2, 3}];
A  =  
        Table[Random[Real, {-3, 3}], {i, 1, hangdim}, {j, 1, 
            hitdim}];
MatrixForm[A]

Now look atMathematica's calculation of PseudoInverse[A]:

MatrixForm[PseudoInverse[A]]

And Mathematica's attempt at a calculation of  ((((A^t) . A))^(-1)) (A^t).

MatrixForm[Inverse[Transpose[A] . A] . Transpose[A]]

The calculation failed. 
The upshot: The formula
               PseudoInverse[A] =  ((((A^t) . A))^(-1)) (A^t).
fails for this matrix.
What gives?

The formula
               PseudoInverse[A] =  ((((A^t) . A))^(-1)) (A^t)
fails for a very good reason: The matrix A is not of full rank.

rank = Length[SingularValues[A][[2]]];
rank  ==  hitdim

T.2) The Spectral Theorem says all symmetric matrices are diagonalizable. You just go with the eigenvectors for your hanger and aligner frames and the eigenvalues for your stretches.

T.2.a.i) When you make your aligner frame the same as your hanger frame, you make a symmetric matrix

You can make a symmetric matrix by:
a) Specifying a perpendicular frame which you use for both your aligner frame and your hanger frame.
b) Going with any stretch factors you like (including positive, zero or negative).
Try it and see what you get in 2D:

Clear[perpframe, alignerframe, hangerframe];
\(\(s = Random[Real, {\(-1.5\), 1.5}];\)
\)
\(\({perpframe[1], perpframe[2]} = 
        N[{{Cos[s], Sin[s]}, {Cos[s + π\/2], 
              Sin[s + π\/2]}}];\)
\)
\(\({alignerframe[1], alignerframe[2]} = {perpframe[1], 
          perpframe[2]};\)
\)
\(\({hangerframe[1], hangerframe[2]} = {perpframe[1], 
          perpframe[2]};\)
\)
\(\({xstretch, ystretch} = {Random[Real, {\(-2\), 2}], 
          Random[Real, {\(-2\), 2}]};\)
\)
aligner = {alignerframe[1], alignerframe[2]};
stretcher = {{xstretch, 0}, {0, ystretch}};
\(\(hanger = Transpose[{hangerframe[1], hangerframe[2]}];\)
\)
A = hanger . (stretcher . aligner);
MatrixForm[A]

A == Transpose[A]

Why did that work? 
Why will this work in any dimension?

Well, 
        A = hanger . stretcher . aligner.
Because the aligner frame is the same as the hanger frame, you are guaranteed that 
        (hanger^t) = aligner
and 
         (aligner^t) = hanger.
And because stretcher is a diagonal matrix,
         (stretcher^t) = stretcher.
So
          (A^t) =    (aligner^t) . (stretcher^t) . (hanger^t)
               = hanger . stretcher . aligner
              = A. 
This tells you that A is symmetric.

T.2.a.ii) The Spectral Theorem tells you that all symmetric matrices can be made this way

How does the Spectral Theorem tell you that all symmetric matrices can be made this way?

\(TextForm\`Go\ with\ any\ symmetric\ matrix\ A\ hitting\ on\ hitdimD . 

        The\ Spectral\ theorem\ gives\ you\ an\ orthonormal\ basis\ of\ hitdimD\ \((perpendicular\ frame\ spanning\ all\ of\ hitdimD)\)\)
           {(Subscript[X, 1]), (Subscript[X, 2]), (Subscript[X, 3]),  
      (Subscript[X, 4]), ...., (Subscript[X, hitdim])} 
\(TextForm\`	\(\(consisting\)\(\ \)\(of\)\(\ \)\(unit\)\(\ \)\(eigenvectors\)\(\ \)\(of\)\(\ \)\(A\)\(\ \)\(with\)\(\ \)\(corresponding\)\(\ \)\(eigenvalues\)\(\ \ \)\)\)
\(TextForm\`	\ \ \ \ \ \ \ \ \ {λ\_1, λ\_2, λ\_3, λ\_4, \(\(...\)\(.\)\), λ\_hitdim}\ \ \)
\(TextForm\`	so\ that\)
\(TextForm\`\ \ \ \ \ \ \ \ \ \ \ A . X\_j\  = \ \(λ\_j\) 
        X\_j\ for\ all\ j' s . 
	When\ you\ put\)
            aligner  =   {(Subscript[X, 1]), (Subscript[X, 2]), (Subscript[X, 3]),  
        (Subscript[X, 4]), ...., (Subscript[X, hitdim])}
\(TextForm\`\ \ \ \ \ \ \ \ \ \ \ \ stretcher\  = \ 
      DiagonalMatrix[{λ\_1, λ\_2, λ\_3, λ\_4, \(\(...\)\(.\)\), λ\_hitdim}\ ]\)
\(TextForm\`\ \ \ \ \ \ \ \ \ \ \ \ \ hanger\  = \ aligner\^t, 

    it' s\ automatic\ that\)
\(TextForm\`\ \ \ \ \ \ \ \ \ \ \ \ \ \((hanger . stretcher . 
            aligner\ )\) . \ X\_j\  = \ \(λ\_j\) 
        X\_j\ \ \ \ \ for\ all\ the\ j' \(\(s . 
Because\)\(\ \)\)\)
\(TextForm\`		\ \ \ \ \ \ \ \ \ \ \ A . 
        X\_j\  = \ \(\(λ\_j\) \(X\_j\)\(\ \ \)\(for\)\(\ \)\(all\)\(\ \)\(the\)\(\ \)\(j'\) \(s\)\(\ \)\)\)
and because  {(Subscript[X, 1]), (Subscript[X, 2]), (Subscript[X, 3]),  (Subscript[X, 4]), ...., 
        (Subscript[X, hitdim])} is a basis for hitdimD,  
    it' s now automatic that
             A  =  
      hanger . stretcher . aligner .

T.2.a.iii) The   Spectral Theorem tells you that all symmetric matrices are diagonalizable

How does the Spectral Theorem tell you that all symmetric matrices are diagonalizable?

\(\(\(Part\ \ ii\)\()\)\)\ above\ said\ that\ when\ you\ go\ with\ any\ symmetric\ matrix\ A\ hitting\ on\ hitdimD, 

      the\ Spectral\ theorem\ gives\ you\ an\ orthonormal\ basis\ of\ hitdimD\ \((perpendicular\ frame\ spanning\ all\ of\ hitdimD)\)\)
           {(Subscript[X, 1]), (Subscript[X, 2]), (Subscript[X, 3]),  
            (Subscript[X, 4]), ...., (Subscript[X, hitdim])} 
\(\(	\)\(consisting\ of\ unit\ eigenvectors\ of\ A\ with\ corresponding\ eigenvalues\)\)  
\(\(	\ \ \ \ \ \ \ \ \ \)\({λ\_1, λ\_2, λ\_3, λ\_4, \(\(...\)\(.\)\), λ\_hitdim}\)\)  
	so that
           \(A . X\_j\  = \ \(λ\_j\) X\_j\ for\ all\ j' 
            s . 
	When\ you\ put\)
            aligner  =   {(Subscript[X, 1]), (Subscript[X, 2]), (Subscript[X, 3]),  
            (Subscript[X, 4]), ...., (Subscript[X, hitdim])}
            \(stretcher\  = \ 
          DiagonalMatrix[{λ\_1, λ\_2, λ\_3, λ\_4, \(\(...\)\(.\)\), λ\_hitdim}\ ]\)
             \(hanger\  = \ aligner\^t, 

        it' s\ automatic\ that\)
        \(A\  = \ 
          hanger . stretcher . aligner . 

              Because\ hanger\ and\ aligner\ are\ both\ based\ on\ the\ same\ frame, \ it' s\ automatic\ that\)
          \(aligner\  = \ hanger\^\(-1\) . 
So\)
           \(A\  = \ 
          hanger . stretcher . hanger\^\(-1\) . 
Now\ put\)
			          SpannerMatrix  =  hanger
              diagonalmatrix  =  stretcher
and read off
\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(A\  = \ 
          SpannerMatrix . diagonal . SpannerMatrix\^\(-1\) . 

              This\ is\ enough\ to\ be\ able\ to\ proclaim\ that\ A\ is\ \(\(diagonalizable\)\(.\)\)\)\)

T.3) Gradient vectors:  gradf[x, 
          y] =  ∇f[x, 
              y]  = {(∂_x) f[x, y], (∂_y) f[x, y]}.
Hessian matrices : (Subscript[H, f])[x, y] =  (({{\(∂\_{x, 2}f[x, y]\), (∂_(x, y))f[x, 
                    y]}, {(∂_(x, y))f[x, y], \(∂\_{y, 2}f[x, 
                    y]\)}})).
Local maximizers, minimizers and saddle points of functions

T.3.a) Background: Gradient vectors and Hessian matrices: 
Using them for  function max-min

Take any function f[x, y].
The gradient of f[x, y] at a point {x, y} is
             gradf[x, 
          y] =  ∇f[x, 
              y]  = {(∂_x) f[x, y], (∂_y) f[x, y]}

Clear[f, gradf, x, y];
gradf[x_, 
        y_] = {(∂_x) f[x, y], (∂_y) f[x, y]}

See the gradient ∇f[x, y] for 
        f[x, y] = (x^2) Sin[y]:

\(\(f[x_, y_]\  = x\^2\ Sin[y];\)
\)
gradf[x, y]

See the gradient ∇f[x, y] for
        f[x, y] = (E^(2 x)) Cos[y]:

Clear[f];
\(\(f[x_, y_]\  = E\^\(2\ x\)\ Cos[y];\)
\)
gradf[x, y]

Get it? 

Take any function f[x,y].
The Hessian matrix (Subscript[H, f])[x, y] of f[x,y] at a point {x,y} is

(Subscript[H, f])[x, y] =  (({{\(∂\_{x, 2}f[x, y]\), (∂_(x, y))f[x, 
                    y]}, {(∂_(x, y))f[x, y], \(∂\_{y, 2}f[x, 
                    y]\)}}))

Clear[f, x, y, H];
	(Subscript[H, f])[x_, y_]= (({{\(∂\_{x, 2}f[x, y]\), (∂_(x, y))f[x, 
                      y]}, {(∂_(x, y))f[x, y], \(∂\_{y, 2}f[x, 
                      y]\)}}))
MatrixForm[
      (Subscript[H, f])[x, y]]

The hessian matrix is guaranteed to be symmetric.

(Subscript[H, f])[x, y]  ==  Transpose[(Subscript[H, f])[x, y]]

See the gradient ∇f[x,y] for and the Hessian (Subscript[H, f])[x, y]  for
              f[x, y] = (x^2)   +  3 x y  + (y^2):

\(\(f[x_, y_]\  = \ 5\ x\^2\ \  + \ 3\ x\ y\  + 4\ y\^2;\)
\)
\(gradf[x, y]
\)
MatrixForm[(Subscript[H, f])[x, y]]

See the gradient ∇f[x,y] for and the Hessian (Subscript[H, f])[x, y]   for
          f[x, y] = (x^2) Sin[y]:

\(\(f[x_, y_]\  = \ x\^2\ Sin[y];\)
\)
\(gradf[x, y]
\)
MatrixForm[(Subscript[H, f])[x, y]]

See the gradient ∇f[x,y] for and the Hessian (Subscript[H, f])[x, y]   for
         f[x, y] = (E^(2 x)) Cos[y]:

Clear[f];
\(\(f[x_, y_]\  = E\^\(2\ x\)\ Cos[y];\)
\)
\(gradf[x, y]
\)
MatrixForm[(Subscript[H, f])[x, y]]

Get it?

The recipe for using the Hessian for function max-min is simple in principle.
You take your function, find {(Subscript[x, 0]), (Subscript[y, 0])} so that
           ∇f[(Subscript[x, 0]), (Subscript[y, 0])] = {0, 0}.
 You are guaranteed that {(Subscript[x, 0]), (Subscript[y, 0])} is 
-> a local minimizer of f[x,y] if both eigenvalues of the Hessian (Subscript[H, f])[(Subscript[x, 0]), (Subscript[y, 0])] are positive.
 
-> a local maximizer of f[x,y] if both eigenvalues of the Hessian (Subscript[H, f])[(Subscript[x, 0]), (Subscript[y, 0])] are negative
 
-> a saddle point of f[x,y] if  (Subscript[H, f])[(Subscript[x, 0]), (Subscript[y, 0])] has one positive and one negative eigenvalue.
 
-> Otherwise, you get no information.

T.3.a.i) Trying it out

Here's a function:

Clear[f, x, y, gradf, H];
f[x_, y_] = 0.4 ((x^2) + 4 Sin[(y^2)])

Analyze some of the points at which the gradient of f[x,y] is {0,0} determining whether they are local maximizers, local minimizers or saddle points.

Here are the gradient and Hessian of f[x,y]:

Clear[f, x, y, gradf, H];
f[x_, y_]  = (1/2) ((x^2) + 4 Sin[(y^2)]);

gradf[x_, 
        y_] = {(∂_x) f[x, y], (∂_y) f[x, y]}
(Subscript[H, f])[x_, y_]=  (({{\(∂\_{x, 2}f[x, y]\), (∂_(x, y))f[
                        x, y]}, {(∂_(x, y))f[x, y], \(∂\_{y, 2}f[
                        x, y]\)}}));

MatrixForm[(Subscript[H, f])[x, y]]

Check out where the gradient is {0,0}:

Solve[gradf[x, y]  ==  0, {x, y}]

This gives rise to several critical \(points . 
Three\) of them are:
{(Subscript[x, 1]), (Subscript[y, 1])} = {0, 0} , {(Subscript[x, 2]), (Subscript[y, 2])} =  {0, (Sqrt[(π/2)])}   and {(Subscript[x, 3]), (Subscript[y, 3])} = {0, -(Sqrt[(π/2)])}.
		Check them out:

{(Subscript[x, 1]), (Subscript[y, 1])}  =  {0, 0};
Eigenvalues[
      (Subscript[H, f])[(Subscript[x, 1]), (Subscript[y, 1])]]

{(Subscript[x, 2]), (Subscript[y, 2])}  =  {0, (Sqrt[(π/2)])};
Eigenvalues[
      (Subscript[H, f])[(Subscript[x, 2]), (Subscript[y, 2])]]

{(Subscript[x, 3]), (Subscript[y, 3])}  =  {0, -(Sqrt[(π/2)])};
Eigenvalues[
      (Subscript[H, f])[(Subscript[x, 3]), (Subscript[y, 3])]]

The Hessian test tells you that  {(Subscript[x, 1]), (Subscript[y, 1])} is a local minimizer  and both {(Subscript[x, 2]), (Subscript[y, 2])} and {(Subscript[x, 3]), (Subscript[y, 3])}  are saddle points.

See a plot of f[x,y] on a circle centered at the local minimizer at {(Subscript[x, 1]), (Subscript[y, 1])}:

Clear[x, y, r, s];
{xbase, ybase}  =  {(Subscript[x, 1]), (Subscript[y, 1])};
x[r_, s_]  =  xbase +  r  Cos[s];
y[r_, s_]  =  ybase +  r  Sin[s];

fplot  =  
          ParametricPlot3D[{x[r, s], y[r, s], f[x[r, s], y[r, s]]}, {r, 0, 
              0.5}, {s, 0, 2 π}, ViewPoint -> CMView, 
            AxesLabel -> {"x", "y", "f[x,y]"}, Boxed -> False, 
            PlotLabel -> "fplot", DisplayFunction -> Identity];

\(Show[fplot, 
        Arrow[{0, 0, 0.5}, Tail -> {xbase, ybase, f[xbase, ybase]}, 
          VectorColor -> Red], 

        DisplayFunction -> $DisplayFunction];\)

Yessiree bob! Just as the Hessian predicted. A local minimizer at:

{(Subscript[x, 1]), (Subscript[y, 1])}

Now check out the saddle point at:

{(Subscript[x, 2]), (Subscript[y, 2])}

Here you go:

Clear[x, y, r, s];
{xbase, ybase}  =  {(Subscript[x, 2]), (Subscript[y, 2])};
x[r_, s_]  =  xbase +  2  r Cos[s];
y[r_, s_]  =  ybase +  r Sin[s];

fplot  =  
          ParametricPlot3D[{x[r, s], y[r, s], f[x[r, s], y[r, s]]}, {r, 0, 
              0.4}, {s, 0, 2 π}, ViewPoint -> CMView, 
            AxesLabel -> {"x", "y", "f[x,y]"}, Boxed -> False, 
            PlotPoints -> {30, 30}, PlotLabel -> "fplot", 
            DisplayFunction -> Identity];

\(Show[fplot, 
        Arrow[{0, 0, 0.5}, Tail -> {xbase, ybase, f[xbase, ybase]}, 
          VectorColor -> Red], 

        DisplayFunction -> $DisplayFunction];\)

Jump on and ride the bronco! Just as the Hessian predicted. A saddle point at:

{(Subscript[x, 2]), (Subscript[y, 2])}

Now check out the saddle point  at:

\(\(\({x\_3, y\_3}\)\(
\)\)\)

Here you go:

Clear[x, y, r, s];
\(\({xbase, ybase}\  = \ {x\_3, y\_3};\)
\)
x[r_, s_]  =  xbase +  2  r Cos[s];
\(\(y[r_, s_]\  = \ ybase + \ r\ Sin[s];\)
\)
\(\(fplot\  = \ 
        ParametricPlot3D[{x[r, s], y[r, s], f[x[r, s], y[r, s]]}, {r, 0, 
            0.3}, {s, 0, 2\ Pi}, ViewPoint -> CMView, 
          AxesLabel -> {"x", "y", "f[x,y]"}, Boxed -> False, 
          PlotLabel -> "fplot", 
          DisplayFunction -> Identity];\)
\)
\(\(Show[fplot, 
        Arrow[{0, 0, 0.2}, Tail -> {xbase, ybase, f[xbase, ybase]}, 
          VectorColor -> Red], 

        DisplayFunction -> $DisplayFunction];\)\)

Yep! Another saddle point at {(Subscript[x, 3]), (Subscript[y, 3])}.
 Just as the Hessian predicted.

T.3.a.ii) Why the Hessian test works

Why does the Hessian test work?

Near any point {a, b}, 
 approxf[x, y]= f[a, b] + ∇f[a, b] . {x  -  a, y  -  b} +{x  -  a,  
                      y  -  b}  . ((((Subscript[H, f])[a, b] . {x  -  a,  
                          y  -  b}))/2)
 mimics the behavior of f[x, y].

Go with a critical point {a,b}. This makes   ∇f[a,b] = {0,0}, and so the gradient term in approxf[x,y] drops out, leaving
        approxf[x, y]= f[a, b]  +{x  -  a,  
                      y  -  b}  . ((((Subscript[H, f])[a, b] . {x  -  a,  
                          y  -  b}))/2)
Now go with {x,y} ≠ {a,b} and put
          {x, y} = s eigvect[1] + t eigvect[2] + {a, b},
so that at least one of s and t is not 0.

This gives
             approxf[x, y]= f[a, b]  +(((s eigvect[1]  +  t eigvect[2],).(Subscript[H, f])[a, b].(s eigvect[1]  +  
                            t eigvect[2],).)/2)

Multiply out on the right to get

approxf[x, y] =           f[a, b]  +(((s eigvect[1]  +  t eigvect[2],).(s eigval[1]  eigvect[1] +  t eigval[2] eigvect[2],).)/2)

This is the same as

approxf[x, y] =     f[a, b]  +((s eigvect[1] .(s eigval[1]  eigvect[1] +  t eigval[2] eigvect[2],).)/2)
                           + (( t eigvect[2],.(s eigval[1]  eigvect[1] +  t eigval[2] eigvect[2],).)/2)

This is the same as

approxf[x, y] =     f[a, b]  +((s eigvect[1] .(s eigval[1]  eigvect[1] ).)/2)
                                          + (( t eigvect[2],.(  t eigval[2] eigvect[2]),)/2).

approxf[x, y] =     f[a, b]  +((s  (s eigval[1]  ).)/2)
                                          + (( t ( t eigval[2] ))/2).

This is the same as
          approxf[x, y] =   f[a, b]  +   (s^2)  eigval[1]  +  
          (t^2) eigval[2].
Now milk this.

Why two positive eigenvalues of (Subscript[H, f])[a, b] reveal that {a, b} is a local minimizer of f[x, y]

If eigval[1] and eigval[2] are both positive, then   (s^2)  eigval[1]  +  (t^2) eigval[
            2]  is also positive (because at least one of s and t is not 0) and so
               approxf[x, y]=  f[a, b] +  (s^2)  eigval[1] + (t^2) eigval[2] 
			                       =   f[a, b] +  positive >  f[a, b] = approxf[a, b]
This means  {a,b} minimizes approxf[x,y]. 
And because approxf[x,y] mimics f[x,y] for {x,y} near {a,b}, this tells you that {a,b}  is a local minimizer of f[x,y].

Why two negative eigenvalues of (Subscript[H, f])[a, b] reveal that {a, b} is a local  maximizer of f[x, y]

If eigval[1] and eigval[2] are both negative, then   (s^2)  eigval[1]  +  (t^2) eigval[
            2]  is also negative (because at least one of s and t is not 0)  and so
              approxf[x, y]=  f[a, b] +  (s^2)  eigval[1] + (t^2) eigval[2] 
=   f[a, b] +  negative <  f[a, b] = approxf[a, b]
This means  {a,b} maximizes approxf[a,b]. 
And because approxf[x,y] mimics f[x,y] for {x,y} near {a,b}, this tells you that {a,b}  is a local maximizer of f[x,y].

Why one positive eigenvalue and one negative eigenvalue of (Subscript[H, f])[a, b] reveal that that the plot of f[x, y] has a saddle point at {a, b, f[a, b]}.

If eigval[1] >0 and eigval[2] < 0  then 
               (s^2)  eigval[1]  +  (t^2) eigval[
            2]  
 is positive for t = 0 and negative for s = 0. (because at least one of s and t is not 0).
When you go with t = 0 and s not zero, you get
                approxf[x, y]=  f[a, b] +  (s^2)  eigval[1] + (t^2) eigval[2]
                     =   f[a, b] +  positive >  f[a, b] = approxf[a, b]

When you go with s = 0 and t not zero, you get
           approxf[x, y]=  f[a, b] +  (s^2)  eigval[1] + (t^2) eigval[2] 
                             = f[xx, yy] +  negative <  f[a, b] = approxf[a, b]

This means that the plot of approxf[x,y] has a saddle point at {a,b}. And because approxf[x,y] mimics f[x,y] for {x,y} near {a,b}, this tells you that the plot of f[x,y] has a saddle point at {a,b}.

This also tells you that 
-> When you leave {a,b} in the direction of eigvect[1], then f[x,y] initially goes up.
-> When you leave {a,b} in the direction of eigvect[2], then f[x,y] initially goes down.

T.3.a.iii) Why 
approxf[x,y] = f[a, b] + ∇f[a, b] . {x  -  a, y  -  b} +((  {x  -  a, 
                y  -  b}  . ((Subscript[H, f])[a, b] . {x  -  a, 
                    y  -  b}))/2)
                             
 mimics f[x,y] near {a,b}

Explain this:
Near any point {a, b}, 
 approxf[x,y] =   \(\(\(f[a, 
          b]\  + \ ∇f[a, b] . {x\  - \ a, 
            y\  - \ b}\  + \(\(\ \ \)\({x\  - \ a, y\  - \ b}\  . \((H\_f[a, b] . {x\  - \ a, y\  - \ b})\)\)\)\/2\)\(
\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\)\)
 mimics the behavior of f[x, y].

Take any fixed point {a, b} and put
       \(approxf[x, y]\  = \ 
          f[a, b]\  + \ ∇f[a, b] . {x\  - \ a, 
                y\  - \ 
                  b}\  + \(\(\ \ \)\({x\  - \ a, y\  - \ b}\  . \((H\_f[a, b] . {x\  - \ a, y\  - \ b})\)\)\)\/2\)
                             
Comparef[a, b]  and approxf[a, b]:

Clear[f, gradf, hessianf, x, y, a, b, xbase, H, approxf, 
        ybase];
gradf[x_, 
            y_] = {(∂_x) f[x, y], (∂_y) f[x, y]};

(Subscript[H, f])[x_, y_]= (Subscript[H, f])[x_, y_]= (({{\(∂\_{x, 2}f[x, y]\), (∂_(x, y))f[
                          x, y]}, {(∂_(x, y))f[x, y], \(∂\_{y, 2}f[
                          x, y]\)}}));
		
approxf[x_, 
            y_]  = 
f[a, 
              b]  +  ∇f[a, b] . {x  -  a, 
                y  -  
                  b}  + ((  {x  -  a, y  -  b}  . ((Subscript[H, f])[a, b] . {x  -  a, y  -  b}))/2) ;
f[a, b] 

approxf[a, b]

\(TextForm\`The\ \(result : 
		f[a, 
            b]\)\ \  = \ \(\(approxf[a, 
              b]\  . 
		Now\ compare\ the\ first\ derivatives\ at\ {a, 
            b}\)\(:\)\)\)

\({∂\_x f[x, y] /. {x -> a, 
          y -> b}, 
	\ \ \ \ \ \ ∂\_x approxf[x, y] /. {x -> a, 
          y -> b}}\)

{(∂_y)f[x, y]/.{x -> a, y -> b},
	      (∂_y)approxf[x, y]/.{x -> a, y -> b}}

The call:
     f[x,y] and approxf[x,y] have the same first derivatives at {a,b}.

\(TextForm\`		\)

Compare second derivatives at with respect to x and y at {a,b}:

{(∂_({x, 2}))f[x, y]/.{x -> a, y -> b},
       (∂_({x, 2}))approxf[x, y]/.{x -> a, y -> b}}

{(∂_({y, 2}))f[x, y]/.{x -> a, y -> b},
         (∂_({y, 2}))approxf[x, y]/.{x -> a, y -> b}}

{(∂_(x, y))f[x, y]/.{x -> a, y -> b},
       (∂_(x,y))approxf[x, y]/.{x -> a, y -> b}}

The upshot: f[x, y] and approxf[x, y] have order of contact 2 at {a, b}.
This tells you that near any point {a, b}, 
      \(\(\(f[a, 
          b]\)\(\ \)\(+\)\(\ \)\(∇f[a, b] . {x\  - \ a, 
            y\  - \ b}\)\(\ \)\(+\)\(\(\(\ \ \)\({x\  - \ a, 
                y\  - \ b}\  . \((H\_f[a, b] . {x\  - \ a, 
                    y\  - \ 
                      b})\)\)\)\/2\)\(
\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\)\)
 mimics the behavior of f[x, y].

T.4) Quadratic forms
         f[x, y]  = a (x^2) + b x y + c (y^2) + d x + e y + g
Ellipses, hyperbolas and parabolas defined by setting a quadratic form equal to a constant

T.4.a.i) A tilted off-set ellipse defined by setting a quadratic form equal to a constant

Here's an example of something folks call a quadratic form:

Clear[f, x, y];
f[x_, y_] = 
      1.1 (x^2) + 1.9 x y + 2.5 (y^2) + 12.4 x + 10.8 y

And look at this plot:

constant = 5;
curveplot = 
        ContourPlot[f[x, y], {x, -15, 5}, {y, -6, 6}, 
          Contours -> {constant}, Axes -> True, 
          AxesLabel -> {"x", "y"}, 
          ContourSmoothing -> Automatic, PlotPoints -> 50, 
          AspectRatio -> Automatic, 
          ColorFunction -> Automatic];

The border of the black region is a plot of the curve consisting of all the points {x,y} for which
             f[x,y] = constant = 5.
Parameterize and plot this ellipse.
Give the perpendicular frame on which it is hung and measure the length of the long axis and the short axis.

Calculate the gradient of f[x,y] and find out where it is {0,0}:

Clear[x, y, gradf, H];
gradf[x_, 
          y_] = {(∂_x) f[x, y], (∂_y) f[x, y]};
\(\(criticals\  = Solve[gradf[x, y]\  == \ 0];\)
\)
{xcritical, ycritical}  =  {x, y} /. criticals[[1]]

Throw a plot of {xcritical, ycritical} into the mix:

embellishedplot  = 
        Show[curveplot, 
          Graphics[{Red, PointSize[0.04], 
              Point[{xcritical, ycritical}]}]];

Now calculate a perpendicular frame of eigenvectors of the Hessian of f[x,y] and throw them into the plot:

Clear[H];
(Subscript[H, f])[x_, y_]=  (({{\(∂\_{x, 2}f[x, y]\), (∂_(x, y))f[
                        x, y]}, {(∂_(x, y))f[x, y], \(∂\_{y, 2}f[
                        x, y]\)}}));

Clear[eigenvector];
{eigenvector[1], eigenvector[2]}  =  
          Eigenvectors[(Subscript[H, f])[x, y]];

scalefactor  =  3;
eigenplot  =  
          Table[Arrow[scalefactor eigenvector[k], Tail -> {0, 0}, 
              VectorColor -> Gold, HeadSize -> 1], {k, 1, 2}];

setup  =  
        Show[embellishedplot , eigenplot];

Put the tails of the eigenvectors at {xcritical,ycritical}:

\(\(adjustedeigenplot\  = \ 
        Table[Arrow[scalefactor\ eigenvector[k], 
            Tail -> {xcritical, ycritical}, VectorColor -> Gold, 
            HeadSize -> 1], {k, 1, 2}];\)
\)
goodsetup  =  
        Show[embellishedplot , adjustedeigenplot];

The natural coordinate system for this set up is the {u,v} coordinate system coming from the two plotted unit eigenvectors. This is the perpendicular frame on which the ellipse is hung.

 To go from uv-coordinates to xy-coordinates, you just use:

Clear[u, v];
{x[u_, v_], y[u_, v_]} = {xcritical, ycritical} + u eigenvector[1] + 
        v eigenvector[2]

Look at this:

Expand[f[x[u, v], y[u, v]]] == constant

This is the same as:

\(Expand[
          f[x[u, v], 
            y[u, v]]]\  - \ \((f[x[u, v], y[u, v]] /. {u -> 0, 
              v -> 0})\)
			 == \ 
      constant\  - \ \((f[x[u, v], y[u, v]] /. {u -> 0, v -> 0})\)\)

Extract the numerical constants:

a = Coefficient[f[x[u, v], y[u, v]], (u^2)]

b = Coefficient[f[x[u, v], y[u, v]], (v^2)]

r = 
      constant - (f[x[u, v], y[u, v]] /. {u -> 0, 
              v -> 0})

Now move in with the u-v parametrization:

Clear[t];
\(	{u[t_], 
        v[t_]}\  = {\(\@r\ Cos[t]\)\/\@a, \(\@r\ Sin[t]\)\/\@b}\)

Now go to the x-y parameterization:

Clear[xx, yy];
{xx[t_], yy[t_]}  = { x[u[t], v[t]], y[u[t], v[t]]}

See it:

\(\(ellipseplot\  = \ 
        ParametricPlot[{xx[t], yy[t]}\ , {t, 0, \ 2\ Pi}, 
          PlotStyle -> {{Thickness[0.01], Red}}, 
          DisplayFunction -> Identity];\)
\)
Show[goodsetup , ellipseplot];

Perfecto.
Make it look pretty:

ustretch  =  ((Sqrt[r])/(Sqrt[a])) ;
vstretch  =  ((Sqrt[r])/(Sqrt[b])) ;
\(\(tail\  = \ {xcritical, ycritical};\)
\)
\(\ \(neweigenplot\  = \ {Arrow[\ ustretch\ eigenvector[1], Tail -> tail, 
            VectorColor -> Gold, 
            HeadSize -> 1], 
		Arrow[\(-\ ustretch\)\ eigenvector[1], 
            Tail -> tail, VectorColor -> Gold, HeadSize -> 1], 
		Arrow[\ 
            vstretch\ eigenvector[2], Tail -> tail, VectorColor -> Gold, 
            HeadSize -> 1], 
		Arrow[\(-\ vstretch\)\ eigenvector[2], 
            Tail -> tail, VectorColor -> Gold, HeadSize -> 1]};\)
\)
Show[embellishedplot , ellipseplot,  neweigenplot ];

The lengths of the short and long axes are:

{2 ((Sqrt[r])/(Sqrt[a]))  , 2 ((Sqrt[r])/(Sqrt[b]))  }

The ellipse is centered at:

{xcritical, ycritical}

T.4.a.ii) A tilted off-set hyperbola defined by setting a quadratic form equal to a constant

Here's another example of something folks call a quadratic form:

Clear[f, x, y];
f[x_, y_] = 
      1.1 (x^2) + 2.9 x y + 1.4 (y^2) + 12.4 x + 10.8 y

And look at this plot:

constant = 12;
\(\(ranger = 30;\)
\)
curveplot = 
        ContourPlot[
          f[x, y], {x, -ranger, ranger}, {y, -ranger, ranger}, 
          Contours -> {constant}, Axes -> True, 
          AxesLabel -> {"x", "y"}, 
          ContourSmoothing -> Automatic, PlotPoints -> 50, 
          AspectRatio -> Automatic, 
          ColorFunction -> Automatic];

The border of the black region is a plot of the curve consisting of all the points {x,y} for which
             f[x,y] = constant = 12.
Parameterize and plot this hyperbola.
Give the perpendicular frame on which it is hung..

Calculate the gradient of f[x,y] and find out where it is {0,0}:

Clear[x, y, gradf, H];
gradf[x_, 
          y_] = {(∂_x) f[x, y], (∂_y) f[x, y]};
\(	\(criticals\  = Solve[gradf[x, y]\  == \ 0];\)
\)
{xcritical, ycritical}  =  {x, y} /. criticals[[1]]

Throw a plot of {xcritical, ycritical} into the mix:

embellishedplot  = 
        Show[curveplot, 
          Graphics[{Red, PointSize[0.04], 
              Point[{xcritical, ycritical}]}]];

Now calculate a perpendicular frame of eigenvectors of the Hessian of f[x,y] and throw them into the plot:

Clear[H];

(Subscript[H, f])[x_, y_]=  (({{\(∂\_{x, 2}f[x, y]\), (∂_(x, y))f[x, 
                      y]}, {(∂_(x, y))f[x, y], \(∂\_{y, 2}f[x, 
                      y]\)}}));
Clear[eigenvector];
{eigenvector[1], eigenvector[2]}  =  
          Eigenvectors[(Subscript[H, f])[x, y]];

scalefactor  =  15;
eigenplot  =  
          Table[Arrow[scalefactor eigenvector[k], Tail -> {0, 0}, 
              VectorColor -> Gold, HeadSize -> 3], {k, 1, 2}];

setup  =  
        Show[embellishedplot , eigenplot];

Put the tails of the eigenvectors at {xcritical,ycritical}:

\(\(adjustedeigenplot\  = \ 
        Table[Arrow[scalefactor\ eigenvector[k], 
            Tail -> {xcritical, ycritical}, VectorColor -> Gold, 
            HeadSize -> 3], {k, 1, 2}];\)
\)
goodsetup  =  
        Show[embellishedplot , adjustedeigenplot];

The natural coordinate system for this set up is the {u,v} coordinate system coming from the two plotted unit eigenvectors. This is the perpendicular frame on which the hyperbola is hung.

 To go from uv-coordinates to xy-coordinates, you just use:

Clear[u, v];
{x[u_, v_], y[u_, v_]} = {xcritical, ycritical} + u eigenvector[1] + 
        v eigenvector[2]

Look at this:

Expand[f[x[u, v], y[u, v]]] == constant

This is the same as:

\(Expand[
          f[x[u, v], 
            y[u, v]]]\  - \ \((f[x[u, v], y[u, v]] /. {u -> 0, 
              v -> 0})\)
			 == \ 
      constant\  - \ \((f[x[u, v], y[u, v]] /. {u -> 0, v -> 0})\)\)

Extract the numerical constants:

a = Coefficient[f[x[u, v], y[u, v]], (u^2)]

b = Coefficient[f[x[u, v], y[u, v]], (v^2)]

r = 
      constant - (f[x[u, v], y[u, v]] /. {u -> 0, 
              v -> 0})

Now move in with the u-v parameterizations:

Clear[u1, v1, u2, v2, t];
	{u1[t_], v1[t_]} ={(((Sqrt[r]) Cosh[t])/((Sqrt[Abs])[a])),(((Sqrt[r]) Sinh[t])/((Sqrt[Abs])[b]))}
	{u2[t_], v2[t_]} ={((-(Sqrt[r]) Cosh[t])/((Sqrt[Abs])[a])),(((Sqrt[r]) Sinh[t])/((Sqrt[Abs])[b]))}

Cosh[t] = (1/2) ((E^t) + (E^(-t)))
Sinh[t] = (1/2) ((E^t) - (E^(-t)))

The reason Cos[t] and Sin[t] are used to parametrize ellipses boils down to:
         x = a (Sqrt[r]) Cos[t]  and  y = b (Sqrt[r]) Sin[t] 
 make 
         ((((x/a)))^2) + ((((y/b)))^2) =  r

Clear[a, b, r, t];
Simplify[(((((a (Sqrt[r]) Cos[t])/a)))^2) + (((((b (Sqrt[r]) Sin[t])/b)))^2) ]

ParametricPlot[{Cos[t], Sin[t]}, {t, 0,  2 Pi}];

The reason Cosh[t] and Sinh[t] are used to parametrize hyperbolas boils down to:
         x = a (Sqrt[r]) Cosh[t]  and  y = b (Sqrt[r]) Sinh[t] 
 make 
         ((((x/a)))^2) - ((((y/b)))^2) =  r:

Clear[a, b, r, t];
Simplify[(((((a (Sqrt[r]) Cosh[t])/a)))^2)-(((((b (Sqrt[r]) Sinh[t])/b)))^2) ]

ParametricPlot[{{Cosh[t], Sinh[t]}, {-Cosh[t], 
            Sinh[t]}}, {t, -2, 2}];

That's why 
-> Sinh[t] is called the hyperbolic Sine of t.
-> Cosh[t] is called the hyperbolic Cosine of t.

Now go to the x-y parameterizations:

Clear[x1, y1, y1, y2];
{x1[t_], y1[t_]}  = { x[u1[t], v1[t]], y[u1[t], v1[t]]}
{x2[t_], y2[t_]}  = { x[u2[t], v2[t]], y[u2[t], v2[t]]}

See it:

{tlow, thigh}  =  {-1, 1};
\(\(branch1plot\  = \ 
        ParametricPlot[{x1[t], y1[t]}\ , {t, tlow, thigh}, 
          PlotStyle -> {{Thickness[0.015], Red}}, 
          DisplayFunction -> Identity];\)
\)
\(\(branch2plot\  = \ 
        ParametricPlot[{x2[t], y2[t]}\ , {t, tlow, thigh}, 
          PlotStyle -> {{Thickness[0.015], Red}}, 
          DisplayFunction -> Identity];\)
\)
Show[goodsetup , branch1plot , branch2plot, 
        PlotRange -> All  ];

See more:

{tlow, thigh}  =  {-2, 2};
\(\(branch1plot\  = \ 
        ParametricPlot[{x1[t], y1[t]}\ , {t, tlow, thigh}, 
          PlotStyle -> {{Thickness[0.01], Red}}, 
          DisplayFunction -> Identity];\)
\)
\(\(branch2plot\  = \ 
        ParametricPlot[{x2[t], y2[t]}\ , {t, tlow, thigh}, 
          PlotStyle -> {{Thickness[0.01], Red}}, 
          DisplayFunction -> Identity];\)
\)
niceplot  =  
        Show[goodsetup , branch1plot , branch2plot, 
          PlotRange -> All  ];

Lookin' good and feelin' good.
Make it look pretty:

a = Coefficient[f[x[u, v], y[u, v]], (u^2)]

b = Coefficient[f[x[u, v], y[u, v]], (v^2)]

r = 
      constant - (f[x[u, v], y[u, v]] /. {u -> 0, 
              v -> 0})

\(\(ustretch\  = \ \@r\/\@Abs[a]\ ;\)
\)
\(\(vstretch\  = \ \@r\/\@Abs[b]\ ;\)
\)
tail  =  {xcritical, ycritical};
\(\ \(neweigenplot\  = \ {Arrow[\ ustretch\ eigenvector[1], Tail -> tail, 
            VectorColor -> Red, 
            HeadSize -> 1], 
		Arrow[\(-\ ustretch\)\ eigenvector[1], 
            Tail -> tail, VectorColor -> Red, HeadSize -> 1], 
		Arrow[\ 
            vstretch\ eigenvector[2], Tail -> tail, VectorColor -> Red, 
            HeadSize -> 1], 
		Arrow[\(-\ vstretch\)\ eigenvector[2], 
            Tail -> tail, VectorColor -> Red, HeadSize -> 1]};\)
\)
Show[embellishedplot,  branch1plot , branch2plot,  
        neweigenplot ];

T.4.a.iii) A tilted off-set parabola defined by setting a quadratic form equal to a constant

Here's yet another example of something folks call a quadratic form:

Clear[f, x, y];
f[x_, y_] = 
      2.0 (x^2) + 6.0 x y + 4.5 (y^2) - 2.0 x - 5.0 y

Look at this plot:

constant = 50;
\(\(ranger = 40;\)
\)
curveplot = 
        ContourPlot[
          f[x, y], {x, -ranger, ranger}, {y, -ranger, ranger}, 
          Contours -> {constant}, Axes -> True, 
          AxesLabel -> {"x", "y"}, 
          ContourSmoothing -> Automatic, PlotPoints -> 50, 
          AspectRatio -> Automatic, 
          ColorFunction -> Automatic];

The border of the black region is a plot of the curve consisting of all the points {x,y} for which
             f[x,y] = constant = 50.
Parameterize, identify and plot this curve. 
Identify the point at the tip on the right.

Copy, paste and edit:
Calculate the gradient of f[x,y] and find out where it is {0,0}:

Clear[x, y, gradf, H];
gradf[x_, 
          y_] = {(∂_x) f[x, y], (∂_y) f[x, y]};
\(\(criticals\  = Solve[gradf[x, y]\  == \ 0];\)
\)
{xcritical, ycritical}  =  {x, y} /. criticals[[1]]

This time, there are no critical points.

In spite of this little set back, calculate a perpendicular frame of eigenvectors of the Hessian of f[x,y] and throw them into the plot:

Clear[H];
\(\(H\_f[x_, 
          y_] = \ {{∂\_{x, 2}f[x, y], ∂\_\(x, y\)f[x, 
                y]}, {∂\_\(x, y\)f[x, y], ∂\_{y, 2}f[x, 
                y]}};\)
\)
Clear[eigenvector];
\(\({eigenvector[1], eigenvector[2]}\  = \ 
        Eigenvectors[H\_f[x, y]];\)
\)
scalefactor  =  20;
\(\(eigenplot\  = \ 
        Table[Arrow[scalefactor\ eigenvector[k], Tail -> {0, 0}, 
            VectorColor -> Red, HeadSize -> 3], {k, 1, 2}];\)
\)
setup  =  Show[curveplot , eigenplot];

This looks promising.
The natural coordinate system for this set up is the {u,v} coordinate system coming from the two plotted unit eigenvectors. This is the perpendicular frame on which the curve is hung.
 To go from uv-coordinates to xy-coordinates, you just use:

Clear[u, v];
{x[u_, v_], y[u_, v_]} = 
      u eigenvector[1] + v eigenvector[2]

Look at this:

uvequation = Expand[f[x[u, v], y[u, v]]] == constant

This exhibits v as a quadratic function of u and reveals that the curve is a parabola!
Solve for v in terms of u:

v[u_]  = ((5.26965 u + 6.5 (u^2) - 50)/1.1094)

This is the u-v parametrization with parameter u.:

Now go to the x-y parameterization:

Clear[xx, yy];
{xx[u_], yy[u_]}  = Expand[{ x[u, v[u]], y[u, v[u]]}]

See it:

{ulow, uhigh}  =  {-2, 2};
\(\(parabolaplot\  = \ 
        ParametricPlot[{xx[u], yy[u]}\ , {u, ulow, uhigh}, 
          PlotStyle -> {{Thickness[0.015], Red}}, 
          DisplayFunction -> Identity];\)
\)
Show[curveplot, parabolaplot];

See more:

{ulow, uhigh}  =  {-5, 4};
\(\(parabolaplot\  = \ 
        ParametricPlot[{xx[u], yy[u]}\ , {u, ulow, uhigh}, 
          PlotStyle -> {{Thickness[0.015], Red}}, 
          DisplayFunction -> Identity];\)
\)
niceplot  =  Show[curveplot, parabolaplot];

Okay.
To try to identify the point at the tip on the right, look at:

v[u]

Set v'[u] = 0 and solve for u:

usol  =  Solve[v'[u]  ==  0, u]

In the u-v coordinate system, the tip is at

{utip, vtip}  =  {u, v[u]} /. usol[[1]]

In the x-y coordinate system, the tip is at

{xtip, ytip}  =  {x[utip, vtip], y[utip, vtip]}

Check:

tipplot  =  
        Graphics[{NavyBlue, PointSize[0.04], 
            Point[{xtip, ytip}]}];
Show[niceplot, tipplot];

Nailed it.
And you're out of here.

T.4.b) Eigenvalues and eigenvectors of the Hessian explain why the plots in part a) turned out the way they did

The facts behind the work above are:
Go with the quadratic form
         \(TextForm\`\(\(	\)\(f[x, y]\  = 
        a\ x\^2 + b\ x\ y + c\ y\^2 + d\ x + e\ y + \(\(g\)\(.\)\)\)\)\)

Find {(Subscript[x, 0]) , (Subscript[y, 0])} so that
               gradf[(Subscript[x, 0]) , (Subscript[y, 0])] = {0, 0}.
Calculate two perpendicular unit eigenvectors, eigvec[1] and eigvec[2], and corresponding eigenvalues,eigval[1] and eigval[2] of the Hessian (Subscript[H, f])[x, y]
When you set 
         f[x,y] = constant
and plot the resulting curve, here's what you can get:

Ellipse:

If eigval[1] and eigval[2] are both positive or are both negative, then you get an ellipse centered at {(Subscript[x, 0]) , (Subscript[y, 0])} framed by  eigvec[1] and eigvec[2] .

Hyperbola:

If  eigval[1] > 0 and eigval[2 ] < 0 or eigval[1] < 0 and eigval[2 ] > 0, then you get a hyperbola centered {(Subscript[x, 0]) , (Subscript[y, 0])} and framed by  eigvec[1] and eigvec[2] .

Line:

If  eigval[1] ≠  0 and eigval[2 ]  = 0 , then you get a line running in the direction of eigvec[2].

Parabola:

If gradf[x , y] = {0, 0} has no solution  then you get a parabola framed by  eigvec[1] and eigvec[2] .

Explain some of these good facts.

Go with 
	f[x, y] =a (x^2)+b x y+c (y^2)+d x+e y+g.
If gradf[x, y] = {0, 0}  has a solution {(Subscript[x, 0]) , (Subscript[y, 0])}, then it turns out that

f[x, y] =  f[(Subscript[x, 0]) , (Subscript[y, 0])] + ∇f[(Subscript[x, 0]) , (Subscript[y, 0])] . {x  -  (Subscript[x, 0]), y  -  (Subscript[y, 0])} +                                                                                              (( {x  -  (Subscript[x, 0]), 
                      y  -  (Subscript[y, 0])}  . ((Subscript[H, f])[(Subscript[x, 0]) , 
                          (Subscript[y, 0])] . {x  -  (Subscript[x, 0]), y  -  (Subscript[y, 0])}))/2).

\(TextForm\`Check\ this\ with\ \(\(Mathematica\)\(.\)\(
\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\)\)

Clear[f, x, y, a, b, c, d , e, g, gradf, 
        H];
f[x_, y_]  = 
        a (x^2) + b x y + c (y^2) + d x + e y + 
          g;
gradf[x_, 
          y_] = {(∂_x) f[x, y], (∂_y) f[x, 
              y]};
(Subscript[H, f])[x_, y_]=  (({{\(∂\_{x, 2}f[x, y]\), (∂_(x, y))f[x, 
                      y]}, {(∂_(x, y))f[x, y], \(∂\_{y, 2}f[x, 
                      y]\)}}));
Simplify[ 
      f[(Subscript[x, 0]) , (Subscript[y, 0])]  +  
        gradf[(Subscript[x, 0]) , (Subscript[y, 0])] . {x  -  (Subscript[x, 0]), 
            y  -  (Subscript[y, 0])}  +   ((  {x  -  (Subscript[x, 0]), y  -  (Subscript[y, 0])}  . ((Subscript[H, f])[(Subscript[x, 0]) , (Subscript[y, 0])] . {x  -  (Subscript[x, 0]), y  -  (Subscript[y, 0])}))/2)]

f[x, y]

So 
f[x, y] =  f[(Subscript[x, 0]) , (Subscript[y, 0])] + ∇f[(Subscript[x, 0]) , (Subscript[y, 0])] . {x  -  (Subscript[x, 0]), y  -  (Subscript[y, 0])} +                                                                                              (( {x  -  (Subscript[x, 0]), 
                      y  -  (Subscript[y, 0])}  . ((Subscript[H, f])[(Subscript[x, 0]) , 
                          (Subscript[y, 0])] . {x  -  (Subscript[x, 0]), y  -  (Subscript[y, 0])}))/2).

And because ∇f[(Subscript[x, 0]) , (Subscript[y, 0])] = {0,0}, this simplifies to

So 
   f[x, y] =  f[(Subscript[x, 0]) , (Subscript[y, 0])] +  (( {x  -  (Subscript[x, 0]), 
                        y  -  
                          (Subscript[y, 0])}  . ((Subscript[H, f])[(Subscript[x, 0]) , (Subscript[y, 0])] . {x  -  (Subscript[x, 0]), 
                            y  -  (Subscript[y, 0])}))/2).

Go with {x,y} ≠ {(Subscript[x, 0]) , (Subscript[y, 0])} and put
          {x, y}=s eigvect[1]+t eigvect[2]+{(Subscript[x, 0]) , (Subscript[y, 0])},
so that at least one of s and t is not 0.

This gives
             f[x, y]= f[(Subscript[x, 0]) , (Subscript[y, 0])]  +(((s eigvect[1]  +  t eigvect[2],).(Subscript[H, f])[(Subscript[x, 0]) , (Subscript[y, 0])].(s eigvect[1]  +  
                            t eigvect[2],).)/2)

Multiply out on the right to get

f[x, y] =   f[(Subscript[x, 0]) , (Subscript[y, 0])]  +(((s eigvect[1]  +  t eigvect[2],).(s eigval[1]  eigvect[1] +  t eigval[2] eigvect[2],).)/2)

This is the same as

f[x, y] =     f[(Subscript[x, 0]) , (Subscript[y, 0])]  +((s eigvect[1] .(s eigval[1]  eigvect[1] +  t eigval[2] eigvect[2],).)/2)
                           + (( t eigvect[2],.(s eigval[1]  eigvect[1] +  t eigval[2] eigvect[2],).)/2)

This is the same as

f[x, y] =     f[(Subscript[x, 0]) , (Subscript[y, 0])]  +((s eigvect[1] .(s eigval[1]  eigvect[1] ).)/2)
                                          + (( t eigvect[2],.(  t eigval[2] eigvect[2]),)/2).

And
      f[x, y] =     f[(Subscript[x, 0]) , (Subscript[y, 0])]  +((s  (s eigval[1]  ).)/2)
                                          + (( t ( t eigval[2] ))/2).

This is the same as
          f[x, 
          y] =   f[(Subscript[x, 0]) , (Subscript[y, 0])]  +   (s^2)  eigval[1]  +  
          (t^2) eigval[2].
Now milk this.

When you set f[x,y] = constant, you get:
            f[(Subscript[x, 0]) , 
            (Subscript[y, 0])]  +   (s^2)  eigval[1]  +  (t^2) eigval[2]  =  
        constant 
This is the same as:
           (s^2)  eigval[1]  +  
          (t^2) eigval[2]  =  constant  -    f[
            (Subscript[x, 0]) , (Subscript[y, 0])]   
In the coordinate system coming from eigvec[1] and eigvect[2] with tails at {(Subscript[x, 0]) , (Subscript[y, 0])}, this gives:
▪Ellipses centered at {s,t} = {0,0}  if eigval[1] and eigval[2] are both positive or both negative.
       When you go back to x-y coordinates, these ellipse are centered at the critical point {(Subscript[x, 0]) , (Subscript[y, 0])} and framed by
       eigvec[1] and eigvec[2].
       
▪ Hyperbolas  centered at {s,t} = {0,0}  if eigval[1] > 0 and eigval[2] < 0 or if eigval[1] < 0 and eigval[2] > 0 
   When you go back to x-y coordinates, these hyperbolas are centered at the critical point {(Subscript[x, 0]) , (Subscript[y, 0])}and framed by
    eigvec[1] and eigvec[2].
    
▪ Lines defined by 
           s = +(Sqrt[((constant - f[(Subscript[x, 0]), (Subscript[x, 0])])/eigval)[1]])   and s = -(Sqrt[((constant - f[(Subscript[x, 0]), (Subscript[x, 0])])/eigval)[1]])
 if eigval[1] ≠ 0 and eigval[2] = 0.When you go back to x-y coordinates, these lines run parallel to eigvec[2].

If there is no solution of ∇f[(Subscript[x, 0]) , (Subscript[y, 0])]  =  {0, 0}, then f[x,y] has no maximum, no minimum and no saddle point.
The result: The plot of f[x,y] is an infinite mountain (or valley) whose cross sections are parabolas. (A more detailed explanation of the parabola case is possible but requires lots of algebra.)