# Optimization, Part A

load(draw)$  set_draw_defaults(line_width=2, grid = [2,2], point_type = filled_circle,
          head_type = 'nofilled, head_angle = 20, head_length = 0.5,
             background_color = light_gray,  draw_realpart=false)$     
 fpprintprec : 5$  ratprint : false$

load ("Econ1.mac");

## Preface

Dowling12A.wxmx uses Maxima to review quadratic form analysis and analysis of
unconstrained optimization of functions of economic variables, following a small part of
Ch. 11 and the first part of Ch. 12 of Introduction to Mathematical Economics (3rd ed), 
by Edward T. Dowling, (Schaum's Outline Series), McGraw-Hill, 2012. This text is a bargain, 
with many complete problems worked out in detail. You should compare Dowling's solutions,
worked out "by hand", with what we do using Maxima here.

A code file Econ1.mac as available in the same section (of Economic Analysis with Maxima),
which defines many Maxima functions used in this worksheet.
Use load ("Econ1.mac");

This worksheet is one of a number of wxMaxima files available in the section
    Economic Analysis with Maxima
on my CSULB webpage.

The main subjects of Ch. 11 and 12  are the use of matrices and determinants in 
economic analysis, such as rank of a matrix, the discriminant, leading principal minors, 
the Jacobian matrix, the Hessian matrix and the bordered Hessian matrix.

Part A includes unconstrained optimization examples. In Part B the use of the bordered 
Hessian matrix for second order conditions in optimization with equality constraints is 
discussed, and also the use of eigenvalues determination of sign definiteness. Also 
discussed in some detail in Part B is the subject of Input-Output Analysis in an open 
economy.

Included in 12A is a discussion of the relations between the Taylor expansion theorems 
and Hessian matrix tests for the nature of an extremum, and a Maxima derivation of the 
sign pattern criteria for leading principal minors of the Hessian matrix (for the cases of 
two and three variables) in analyzing unconstrained critical points, which are not part of 
Dowling's text.

We have slightly changed some of the symbols used by Dowling in particular problems.

Quoting Math 2640, Introduction to Optimization (Leed's Univ.)
(http://webprod3.leeds.ac.uk/catalogue/dynmodules.asp?Y=202021&M=MATH-2640)

"Optimisation ''the quest for the best'' plays a major role in financial and economic theory, 
eg., in maximising a company's profits or minimising its production costs. How to achieve 
such optimality is the concern of this course, which develops the theory and practice of 
maximising or minimising a function of many variables, either with or without constraints.
This course lays a solid foundation for progression onto more advanced topics, such as 
dynamic optimisation, which are central to the understanding of realistic economic and 
financial scenarios."

Edwin L. (Ted) Woollett
https://home.csulb.edu/~woollett/
June 13, 2022

## References

Notes on a High School course on optimization methods with a large variety of
approaches to the problem.

https://web.stanford.edu/group/sisl/k12/optimization/#!index.md
-----------------------------------------------------
Course notes, problem sets, and Exams for Math 2640, Introduction to Optimization, 
Leeds Univ., Chris Jones, 2005

http://www1.maths.leeds.ac.uk/~cajones/math2640/MATH2640.html
--------------------------------------------------
Fundamental Methods of Mathematical Economics, Alpha C. Chiang and 
Keven Wainwright, 4th ed., 2005, McGraw-Hill
----------------------------------
Wainwright's course notes at British Columbia Institute of Technology, Burnaby, 
British Columbia, Canada
http://faculty.bcitbusiness.ca/kevinw/chiang/ChapterLectureNotes.htm
----------------------------------------------
Wainwright's 2007 Econ 331 course at Simon Fraser Univ.
http://www.sfu.ca/~wainwrig/Econ331/331.htm
-------------------------------------------------------
Martin J. Gander, Math. Dept.,  Univ. of Geneva,  Ch. 4, Optimization
https://www.unige.ch/~gander/teaching/polycopie.pdf
https://www.unige.ch/~gander/
------------------------------------------------------
Math 2070, Univ. of Sydney Optimization notes

https://www.maths.usyd.edu.au/u/UG/IM/MATH2070/r/NLoptWC.pdf
----------------------------------------------------
Wolfram summary of Mathematics methods for optimization.
https://reference.wolfram.com/language/tutorial/ConstrainedOptimizationIntroduction.html

## rank (amatrix), Linear Independence  [11.1]

If the determinant of a matrix equals zero, the determinant is said to "vanish" and the
matrix is termed "singular".

A "singular matrix" is one in which there exists linear dependence between at least two
rows or columns. 

If the determinant of a matrix is not equal to zero, the matrix is said to be "nonsingular", 
and all its rows and columns are linearly independent. 

If linear dependence exists in a system of equations, the system as a whole will have
an infinite number of possible solutions, making a unique solution impossible.

The rank of a matrix is defined as the maximum number of linearly independent rows
or columns in the matrix. Given a square matrix of order n,
    if rank (A) = n, then A is nonsingular, |A| # 0,  and there is no linear dependence.
    if rank (A) < n, then A is singular, |A| = 0, and there is linear dependence.

The Maxima function rank (M) computes the rank of the square matrix M.
The Maxima function determinant(M) calculates the determinant of M.

### Example 1

A : matrix ( [6, 4], [7, 9] );
B : matrix ( [4, 6], [6, 9] );

The individual elements of a matrix M can be found using M[i,j], in which i and j
are positive integers.

A[1,2];

The elements of row 1 of matrix A is returned as a list, using A[1].

A[1];

list_matrix_entries (M) returns a list of all the elements of the matrix M.

list_matrix_entries(A);

The determinant of a 2 x 2 matrix is called a "second order determinant."
The determinant is only defined for a square matrix.

determinant (A);
determinant (B);

Since |A| # 0, the matrix A is nonsingular and there is no linear dependence between
any of its rows and columns.

rank (A);

The rank of A is 2 which is equal to its dimensions, indicating a nonsingular matrix.

Since |B| = 0, the matrix B is singular and linear dependence exists between its rows
and columns. In this example, row 2 = 3/2 times row1, and col 2 = 3/2 times col 1.

rank (B);

The rank of B is 1 which is less than the dimensions of B (n = 2), and since rank(B) = 1,
there is only one linearly independent row and one linearly independent column in B.

## The Minors of a Matrix   Minor (M, i, j)   [11.3]

We bind the symbol A now to a 3 x 3 matrix with elements a[i,j] using the Maxima
function genmatrix. The a[i,j] elements are components of a Maxima "hash array."
The elements a[i,j] are not (initially) bound to any specific values.

A : genmatrix (a, 3, 3);

The (display2d = false) appearance of this result can be revealed using the Maxima
function grind.

grind(%)$

The Maxima function minor (M, i, j) returns a submatrix of M gotten by removing from M 
row i and column j. In conventional mathematics terminology, a "minor of a matrix" is not 
a submatrix but rather a determinant of the submatrix. Dowling p. 226 uses the term 
"subdeterminant" for a minor of a matrix.

Here we get the matrix which results if we remove row 1 and column 1 from A.

minor (A, 1,1);

grind(%)$

The Maxima function Minor (M, i, j), defined in Econ1.mac,  conforms with conventional 
math terminology (and Dowling), producing a scalar by taking the determinant of a 
submatrix:

fundef (Minor);

Here we find the determinant of the submatrix defined by deleting row 1 and 
column 1 of the 3 x 3 matrix A, using our Maxima function Minor.

Minor (A, 1, 1);

grind(%)$

Minor (A, 3, 3);

## Leading Principal Minors of a Matrix LPM (M, j)  [12.2]

For an n x n square matrix there are n "leading principal minors".

A basic minor of a matrix is the determinant of a square matrix that is of maximal size
with nonzero value. For an n × n nonsingular square matrix, there are n leading 
principal minors.

For the matrix A, a square 3 x 3 matrix, there are 3 "leading principal minors" which
we can obtain using our Maxima function LPM (amatrix, num), defined in Econ1.mac.

The leading principal minor A1 is simply the element A[1,1], produced by deleting 
every row but the first, and deleting every column but the first.

A[1,1];
LPM (A, 1);

The leading principal minor A2 is the determinant of the submatrix produced by 
deleting all but the first 2 rows and columns of A.

LPM (A, 2);

The leading principal minor A3 is the determinant of the matrix produced by deleting all
but the first 3 rows and columns of A, which, since A is a 3 x 3 matrix, is simply the
whole matrix A, since A is defined as a 3 x 3 matrix.

LPM (A, 3);

determinant (A);

## jacobian (funcList, varList), Linear Indepence  [12.1]

Maxima has the function
        jacobian (funcList, varList)
which computes the Jacobian matrix, which can be used to test for functional independence,
 both linear and nonlinear. A Jacobian matrix is composed of all the first-order partial 
derivatives of a system of equations, arranged in an ordered sequence. 
    
 jacobian (funcList, varList) returns the Jacobian matrix of the list of functions funcList 
with respect to the list of variables varList. The (i, j)-th element of the Jacobian matrix is 
diff ( funcList[i],  varList[j] ).

Here is a symbolic example, using Maxima, for  jacobian ( [f1, f2, f3], [x1, x2, x3] ).

depends ( [f1, f2, f3], [x1, x2, x3] );

J : jacobian( [f1, f2, f3], [x1, x2, x3]);

grind(J)$

The first row elements are the first derivatives of f1 with respect to all three variables.
The first column elements are the first derivatives of all the functions with respect to x1.

If we are given a set of three equations
     f1 = 0, f2 = 0, f3 = 0,
then if the determinant of the jacobian J (just defined) is NOT equal to zero, the three 
equations ARE functionally independent. 

If the determinant of the jacobian J IS equal to zero, then the equations are NOT 
independent.

If the list of equations are actually a list of the first derivatives of an expression with
respect to the independent variables, the determinant of the Jacobian matrix is the same 
as the determinant of the Hessian matrix (vide infra), which computes all second order 
derivatives of the given expression, and the latter determinant is, in turn, the same as 
the "n'th leading principal minor" of that Hessian n x n matrix  (where the number of 
independent variables is n). We will approach these concepts slowly.

### Example 1

Given the pair of equations
  f1 = 5 x1 + 3 x2 = 0,
  f2 = 25 x1^2 + 30 x1 x2  + 9 x2^2 = 0,
we first set up the jacobian matrix of first derivatives of f1 and f2.

f1 and f2 are defined as Maxima expressions (rather than Maxima functions).
The implied equations are f1 = 0 and f2 = 0.

f1 : 5*x1  + 3*x2;
f2 : 25*x1^2 + 30*x1*x2 + 9*x2^2;
J : jacobian ([f1, f2], [x1, x2]);

f1 : 5*x1  + 3*x2;
f2 : 25*x1^2 + 30*x1*x2 + 9*x2^2;
J : jacobian ([f1, f2], [x1, x2]);

grind(J)$

Then we use the Maxima functions determinant and expand to evaluate J, the
determinant of the jacobian matrix, (Dowling refers to this as the "Jacobian 
determinant") symbolically.

JD : determinant (J), expand;

Since the determinant of the jacobian evaluates to zero, the two equations f1 = 0, and
f2 = 0, are not functionally independent. The functional dependence is easily seen if
we examine f1^2.

f1^2, expand;
f2;

## Quadratic Forms and the Discriminant  [12.3]

Determinants (in the form of leading principal minors) may be used to test for positive 
or negative definiteness of any quadratic form. The determinant |D| of a quadratic form 
is called a "discriminant".

The matrix D referred to here is a symmetric matrix such that a quadratic form q in an 
arbitrary number of variables can be written as 
         q = transpose(u) . (D u)
in which u is written as a matrix column vector.

The Maxima function cvec ( alist ), defined in Econ1.mac, creates a column matrix whose
elements are the elements of the list alist.

u : cvec ([x, y]);

grind (u)$

A quadratic form q involving two variables (x,y) might be given as
  q = a x^2 + b x y + c y^2,
and we can find a symmetric matrix D by placing the coefficients of the square terms on
the principal diagonal and dividing the coefficients of the nonsquared term equally 
between the off-diagonal positions:

D : matrix ( [a, b/2], [b/2, c]);

transpose (u);

transpose (u) . D . u, expand;

We then evaluate the leading principal minors using our Maxima function LPM(D,j). 
Since D is a 2 x 2 matrix, there are two leading principal minors.

D1 : LPM (D,1);

D2 : LPM (D,2);

If D1 > 0 and D2 > 0, then q is positive definite and q is positive for all values of the 
variables (x,y) as long as x and y are not both zero.

If D1 < 0 and D2 > 0, then q is negative definite and q is negative for all values of the 
variables (x,y) as long as x and y are not both zero.

Note that D2 is just the determinant of D. If D2 is not greater than zero, q is not 
"sign-definite" and q may assume both positive and negative values.

### Two Variable Quadratic Form LPM Test

We can derive these conditions for the definiteness of a two variable quadratic form, 
starting with our general two variable quadratic form expression 
    q =   a x^2 + b x y + c y^2, 
by adding and subtracting a term proportional to y^2 with a coefficient chosen to allow 
us to write q in the form
     q = a (x + d y)^2  + e y^2.
We need to find the values of d and e which makes these alternative expressions
equal.

Let ex0 be the first form, ex1 be the second form, and expr being the difference.
In the following, we want expr to equal zero.

ex0 : a*x^2 + b*x*y + c*y^2;
ex1 : a*(x + d*y)^2 + e*y^2;
expr : ex0 - ex1, expand;

facsum (expr, y^2, x*y);

To arrive at expr = 0, we need the coefficient of y^2 to equal zero, and separately, we want the
coefficient of x*y to equal zero, which gives us two equations in two unknowns.

ey2 : ratcoeff (expr, y^2);

exy : ratcoef (expr, x*y);

soln : solve ([ey2, exy], [d, e] );

soln : soln[1];

e : at( e, soln);

Writing the two variable quadratic form q in the form of ex1

ex1;

q is definitely positive if x and y are not both equal to zero, and if a > 0 and if e > 0.
But a = D1, and e is D2/D1, as we show here:

[D1, D2, expand(D2/D1), expand (e)];

So the general two variable quadratic form q can be written as
      q = D1 (x + d y)^2 + (D2/D1) y^2.

Since the symbols e, D1, and D2 are already bound to expressions, let's kill the binding 
of all three first. 

kill(e, D1, D2)$
subst ([a = D1, e = D2/D1], ex1);

In this form, it is clear that (with x and y not both zero),
     q is positive definite (PD) if D1 > 0 and D2 > 0,
     q is negative definite (ND) if D1 < 0 and D2 > 0.

### Example 3  (2 variable quadratic form)

Given the quadratic form z with numerical coefficients
        z = 2 x^2 + 5 x y + 8 y^2,
form the symmetric matrix D as described above and test for sign-definiteness of
the given quadratic form.

D : matrix ( [2, 5/2], [5/2, 8] );

The leading principal minor D1 is the determinant of the matrix produced by deleting all
but the first row and first column, which determinant is just the matrix element D[1,1] = 2. 

determinant( matrix ([2]));

Using our Maxima function LPM (D, j) for the j'th leading principal minor of the matrix
D, we get:

LPM (D, 1);

The leading principal minor D2 is the determinant of the matrix produced by deleting all
but the first two rows and the first two columns of the matrix D, which will produce the
determinant of D

determinant (D);

LPM (D, 2);

%, numer;

Since both D1 and D2 are greater than zero, z is positive definite and z is positive for
all values of the variables (x,y), excepting the case in which both variables are zero.

### Qtest (amatrix)

The Maxima function Qtest, defined in Econ1.mac, tests a symmetric numerical matrix , and
prints a list of lists: [ ["LPM1", LPM1], ["LPM2", LPM2],...,["LPMn", LPMn] ]. with the 
values of the leading principal minors of the given matrix.

Qtest (D);

### Three Variable Quadratic Form LPM Test

Start with general three variable quadratic form 
     q = a x^2 + b x y + c x z + d y^2 + e y z + f z^2,
which can be written as a matrix equation with D a symmetric 3 x 3 matrix of coefficients
and u a matrix column vector with elements (x,y,z), such that
  q = transpose(u) . D . u

kill(a,b,c,d,e,f)$
D : matrix ( [a, b/2, c/2],  [b/2, d, e/2], [c/2, e/2, f] );

u : matrix ( [x], [y], [z] );

transpose (u) . D . u, expand;

which equals our starting expression for a third order quadratic q.

We now calculate the three leading principal minors of the matrix D.

D1 : LPM (D, 1);

D2 : LPM (D, 2);

D3 : LPM (D, 3);

We now want to prove  that (with x, y and z not all zero)
     q is positive definite (PD) if D1 > 0 and D2 > 0 and D3 > 0,
     q is negative definite (ND) if D1 < 0 and D2 > 0, and D3 < 0.

We look for g, h, k, l, and m (5 unknowns) such that we can write q as a sum of
squares in the form (with a, b, c, d, e, f taken as given):

  q = a (x + g y + h z)^2  + k (y + l z)^2 + m z^2

If a, k, and m are all positive, q is positive for any values of (x,y,z) (not all zero),
    and q is positive definite.
If a, k, and m are all negative, q is negative for any value of (x,y,z) (not all zero),
    and q is negative definite.

Let ex0 be the general starting form of q, ex1 be the alternative expression of
q, and expr be the difference between ex0 and ex1.
 
We require expr to be equal to zero in order to determine g, h, k, l, and m.

kill (g, h, k, l, m)$
ex0 : a*x^2 + b*x*y + c*x*z + d*y^2 + e*y*z + f*z^2;
ex1 : a*(x + g*y + h*z)^2 + k*(y + l*z)^2 + m*z^2;
expr : ex0 - ex1, expand;

We want ex0 and ex1 to be equal, hence expr equal to zero. Hence the coefficient of
z^2 to equal zero,  the coefficient of x*y to be zero, etc.

ez2 : coeff (expr,z^2);

exy : ratcoef (expr, x*y);

exz : ratcoef (expr, x*z);

ey2 : coeff (expr, y^2);

eyz : ratcoef (expr, y*z);

We then ask solve to take the five equations, ez2 = 0, exy = 0, etc and
come up with symbolic solutions for the five unknowns.

solns : solve ([ez2, exy, exz, ey2, eyz], [g, h, k, l, m] );

soln : solns[1];

We bind the symbol k to the result returned in soln.

k : at (k, soln);

We want to show that k is the same as D2/D1.

[D1, D2, expand (D2/D1), expand (k) ];

We see that k = D2/D1.

Next bind the symbol m to the result returned in soln.

m : at (m, soln);

We want to show that m is the same as D3/D2.

[ expand (D3/D2), expand (m) ];

is (equal (expand (D3/D2),  expand (m) ) );

So the general three variable quadratic can be written as
 q = D1 (x + g y + h z)^2 + (D2/D1) (y + l z)^2 + (D3/D2) z^2

Since the symbols k, m, D1, D2, and D3 are already bound to expressions, let's kill 
the binding of all five first.

ex1;

kill (k, m, D1, D2, D3)$
subst ([a = D1, k = D2/D1, m = D3/D2], ex1);

In this form, it is clear that (with x, y and z not all zero)
     q is positive definite (PD) if D1 > 0 and D2 > 0 and D3 > 0,
     q is negative definite (ND) if D1 < 0 and D2 > 0, and D3 < 0. 
  QED.

### Problem 12.6 (3 variable quadratic form)

Given the quadratic form in three variables (x1, x2, x3)
    y = 5 x1^2 - 6 x1 x2 + 3 x2^2 - 2 x2 x3 + 8 x3^2 - 3 x1 x3,
check for sign definiteness.

The coefficients of the squared terms continue to go on the principal diagonal, while the
  coefficient of x1 x2 is equally divided between D[1,2] and D[2,1], etc.

D : matrix ( [5, -3, -3/2], [-3, 3, -1], [-3/2, -1, 8] );

D1 : LPM (D, 1);

D2 : LPM (D, 2);

D3 : LPM (D, 3);

%, numer;

D1, D2, and D3 are all positive and thus the quadratic form y is positive definite and y
is positive for all values of (x1, x2, x3) such that not all three are simultaneously zero.

Note that D3 is the same as the determinant of D.

Qtest (D);

### Problem 12.7 (a) (3 variable quadratic form)

Use discriminants to determine the sign definiteness of the following function:

y = - 2 x1^2 + 4 x1 x2 - 5 x2^2 + 2 x2 x3 - 3 x3^2 + 2 x1 x3.

D : matrix ( [-2, 2, 1], [2, -5, 1], [1, 1, -3] );

D1 : LPM (D, 1);

D2 : LPM (D, 2);

D3 : LPM (D, 3);

Since D1 < 0, D2 > 0, and D3 > 0, y is a negative definite quadratic form and 
y < 0 for all values of x1, x2, x3 as long as they are not all zero simultaneously.

Qtest (D);

## Taylor's Theorem and the Hessian Matrix

In 1 dimension, Taylor's theorem, expanding about the point x, and with dx a vanishingly
small scalar of either sign, is

  f (x + dx) = f(x) + dx fx + (1/2!) dx^2 fxx  + ...

In n dimensions, Taylor's theorem is, with x standing for variables
(x1, x2, ..., xn) and dx standing for (dx1, dx2,....,dxn), and in the following
the third term assumes dx is a matrix column vector:

f(x + dx) = f(x) + dx . grad(f) + (1/2!) transpose(dx) . (H dx) + ...,
where H is the n x n Hessian matrix with components
   H[i,j] = diff (f, xi,1,xj,1) if i # j,
   H[i,j] = diff (f, xi, 2) if i = j.

The general Taylor series formula is derived by considering f(x + λ dx) as a function of λ,
and expanding about λ = 0 using the Taylor expansion formula for a function of one
variable, λ, together with the chain rule. Then putting λ = 1 gives the above expansion
theorem in n dimensions.

For n = 2, we have the two variable case of the Taylor expansion theorem, expanding
about the point (x,y), and with h and k vanishingly small scalars of either sign,

f(x + h, y + k) = f(x, y) + h ∂f/∂x + k ∂f/∂y + 
                        (1/2) ( h^2 ∂²f/∂x² + 2 h k ∂²f/∂x∂y + k^2 ∂²f/∂y² ) + ...

We can write the quadratic term in terms of the Hessian matrix, making use of 
  Young's theorem: ∂²f/∂x∂y = ∂²f/∂y∂x. In the following we use the Maxima function
   hessian (f, [x1, x2]) for a function of two variables. In the following, f is not bound to
   any Maxima expression or Maxima function, so we simply tell Maxima that there
   is some theoretical dependence of f on x and y.

[f, x, y];

depends (f, [x, y]);

That allows us to construct a symbolic 2 x 2 Hessian matrix, which we call H here.

H : hessian (f, [x, y]);

The 2 x 2 Hessian matrix is composed of the second order direct partials (like fxx) on 
the principal diagonal, and the second order cross partials (like fxy) off the principal 
diagonal. Maxima automatically uses the general result that ∂²f/∂x∂y = ∂²f/∂y∂x, and
H is then obviously symmetric.

Let dX be a column matrix with elements (h,k):

[h, k];

dX : matrix ([h], [k] );

The transpose of dX is a row vector with elements (h, k).

transpose (dX);

We can multiply conformable Maxima matrices using spaces and a single period  (.)

transpose (dX) . dX;

transpose (dX) . H . dX;

expand(%);

which is (except for the overall multiplier (1/2)) the quadratic term in the Taylor
expansion of f(x+h, y+k).

We can also show the third order symbolic Hessian matrix form.

z;

depends (f, [x,y,z])$
hessian (f, [x,y,z]);

Notice that the Maxima function hessian automatically uses Young's theorem: 
   ∂²f/∂y∂x = ∂²f/∂x∂y.

## Unconstrained Optimization

If there are no constraints, we have the case "free optimization". 
Let f(x1, x2, ...., xn) be the "objective function" of n variables. 

### optimum (f, varL),  optimumAll (f, varL)

The quickest route to finding unconstrained optimum values is to use optimum(f,varL) 
defined in our software file Econ1.mac. optimum screens out values returned by solve that
include imaginary parts, and given that solve returns real values, screens out real 
negative values. optimum is really meant for Economics problems in which the search 
variables are real and non-negative.

If you want to accept *all* solutions returned by solve, then use optimumAll, with
the same syntax otherwise. See our examples.

In more detail, optimum (f, [x1,x2,....,xn]);  or optimumAll (f, [x1,x2,...,xn]);

### Do it "by hand"...

Otherwise, do it "by hand" (with some help from our other functions).

1. Find the point(s) for which grad(f) = 0. A list of the expressions diff(f, xj) can be
found using gradf : jacobian ([f], [x1, x2, ..., xn])[1];

If the derivatives involve simple polynomials, you can then use Maxima's solve function
as in:
  solns : solve (gradf, [x1, x21, x3,..., xn]).
(If solve has difficulties, or you want to check on solve, you can try to_poly_solve.)

If solve is successful in finding some or all of the solutions, you will get a list of a list, 
as in [ [x1 = 2, x2 = -1/2,...]]  for a single critical point, which you can extract using
   cp : solns[1]  --> [x1 = 2, x2 = -1/2,...].

If solve returns a list of , say, two critical points where the first derivatives of f vanish, such
as    [ [x1 = 2, x2 = -1/2,...], [x1 = 5, x2 = -4,....] ], you could extract these two critical
points using  [cp1, cp2] : solns;  and then cp1 is bound to the list [x1 = 2, x2 = -1/2,...], 
and cp2 is bound to the list [x1 = 5, x2 = -4,....].

2. You can then find the value of the objective function f at each of these points
using
    f1 : at (f, cp1),
    f2 : at (f, cp2).
That will quickly give you a clue as to which is the "biggest" and which is the "smallest"
value of f.

3. To see if a candidate critical point cp satisfies the second order conditions for a relative
maximum or minimum, 
   a.) Calculate the Hessian matrix of f, using the Maxima function hessian (expr, varList),
for example,
H : hessian (f, [x1, x2, ..., xn])  -->  n x n square matrix of second derivatives of f.
   
  b.)  H could turn out to be a purely numerical matrix. In that case, proceed with step c.
Otherwise, we recommend that you turn H into a purely numerical matrix which will depend
on the numerical values of the critical point being investigated.

For example
     Hcp1 : at (H, cp1);

  c.)   Calculate the "leading principal minors" of the purely numerical Hessian matrix Hcp1.
You can use the Maxima function LPM (M, k) for k = 1, 2,..., n (defined in Econ1.mac).
For example 
H1 : LPM (Hcp1, 1);
H2 : LPM (Hcp1, 2);
.....
Hn : LPM (Hcp1,n);

The leading principal minors are the determinants of submatrices of Hcp1. 

  d.) A sufficient condition for a relative MINIMUM is that all the determinants 
  Hk = LPM (Hcp1, k) > 0. Then d²f is a positive definite quadratic form (PD).

  e.) A sufficient condition for a relative MAXIMUM is that (-1)^k*Hk > 0 for all k = 1,2,...,n.
In other words, H1 < 0, H2 > 0, H3 < 0, ...., (-1)^n Hn > 0, and the signs of the 
determinants alternate. Then d²f is negative definite quadratic form (ND).

### Sufficient Conditions vs. Necessary Conditions

We quote from Ch. 1 of Optimization Methods in Economics by John Baxley, retired 
professor in the Mathematics Dept. at Wake Forest Univ.

    https://users.wfu.edu/~baxley/m254book.pdf

Max-min problems play a central role in every calculus course. Finding relative (local) 
maxima and minima using the derivative and applying the first or second derivative test
is the name of the game in curve-sketching as well as the “applied” problems in the 
calculus books.

The student who comes to economics from such calculus courses often feels betrayed. 
Slowly it becomes evident that economists do not spend their time finding maxima and 
minima. In fact, quite the opposite is true. Unlike the typical math problem where one 
“finds the maximum”, the economist assumes that the economic agent (firm, consumer, 
etc.) is instinctively maximizing. The fundamental assumption is that somehow such 
economic agents have a built in computer or natural instinct which leads them to 
maximizing behavior. The central question for the economist is not: find the maximum, 
but: how will the agent adjust maximizing behavior if some variable which he cannot control 
undergoes a change. 

For example, how will the quantity of snack crackers sold in the marketplace change if the 
price of a related good like Coca-cola rises? This question assumes that consumers are 
maximizing their utility and as they face higher prices for Coke, they will make adjustments 
in their expenditures which may effect the amount of snack crackers which are sold.

One of the interesting sidelights of this state of affairs is that economists deeply wish that
second order conditions were necessary (rather than sufficient) for an optimum. You will 
recall the second derivative test:

If c is a critical point of a function f(x) (i.e. f'(c) = 0), then f''(c) < 0 is a sufficient condition 
for c to be a maximum. However, this condition is not necessary, for f(x) = 9 − x^4 has a 
maximum at the critical point x = 0 but f''(0) = 0. 

If we are math students and are on the prowl for maxima, the second derivative test can 
be used to determine if c is a maximum, but if we are economists we want the thought to 
flow the other way: if we know that c gives a maximum, we would like to conclude that 
f''(c) < 0. Unfortunately, we can only conclude that f''(c) ≤ 0. 

The possibility that f''(c) = 0 is often disastrous for economic analysis, because in the 
analysis this value occurs in a denominator and leads to division by zero. Wishful thinking 
has led many economists to argue that this disaster is somehow very unlikely and can be 
safely ignored; some economists have actually referred to the disaster as “pathological” in 
nature. The example f(x) = 9−x^4  does not look very pathological!

As we proceed, you will find that we will often have to assume that second order conditions
known only to be sufficient actually hold at the maximum. We don’t really have a viable
choice. Either we make this assumption and draw an interesting economics conclusion, or
we don’t make the assumption and no conclusion can be drawn.

### Example 2: Extremum of f(x, y)  [12.2]

Let z be a Maxima expression corresponding to:  z = 3 x^2 - x y + 2 y^2 - 4 x - 7 y + 12.
Find any critical points and use the Hessian matrix test to examine the nature of such.

I will use f instead of z.

f : 3*x^2 - x*y + 2*y^2 - 4*x - 7*y + 12;

### optimum(f,varL)

The quickest route is optimum (f,[x,y]).

optimum (f, [x, y] );

optimum has found one (real and non-negative) critical point (x=1,y=2) at
which f has a local minimum. For this kind of problem, optimum calls solve and CPtest.
The solutions returned by solve are screened to avoid solutions which include
imaginary numbers, and given real numbers, avoids non-negative values for (x,y).

optimum defines a global list cp.

cp;

cp[1];

at (f, cp[1]);

### Do it "by hand"...

Let's now take the "long way home" and do things more or less "by hand".
We first let gradf be a list of the first derivatives of f wrt x and y. The fastest way
is to use the Maxima function jacobian, as shown here:

gradf : jacobian ( [f], [x, y])[1];

solns : solve (gradf, [x, y]);

Maxima's solve function has found one critical point. Let's call the "replacement rules"
cp1.

cp1 : solns[1];

What is the value of the Maxima expression f at the critical point?

You can either use

subst (cp1, f);

or 

at (f, cp1);

or the interactive work form (not recommended by the Maxima experts for deep reasons):

f, cp1;

Next we check the second order conditions (SOC) for a relative optimum, given a
critical point which satisfies the first order condition (FOC) that the first derivatives
(wrt x and y) of f simultaneously vanish.

Define the Hessian matrix H, using the Maxima function hessian (expr, varList).

H : hessian (f, [x, y] );

The Hessian test first requires that the Hessian matrix be evaluated for x and y at a
critical point (where the first derivatives vanish). Since for this example, the Hessian
matrix is already purely a set of numbers, we can skip that step.

For any two variable expression or function f(x,y) the Hessian test is that we 
have a relative minimum if the two "leading principal minors" H1 and H2 are both positive. 

We can use our Maxima function LPM (amatrix, num) to return these values,  which 
need to be reduced to numerical values for our test.

We have a relative minimum if the first leading principal minor H1 > 0 and the second
leading principal minor H2 > 0 as well.

We have a relative maximum if the first leading principal minor H1 < 0 and the
second leading principal minor H2 > 0.

We call H1 the first leading principal minor; it is the minor gotten by taking the determinant 
of the matrix left after deleting all rows except the first row and deleting all columns except 
the first column, which leaves a matrix with only one element H[1,1] = fxx = 6. The 
determinant of such a one element matrix is simply the value of H[1,1] = 6.

H1 : LPM (H, 1);

We call H2 the second leading principal minor; the determinant of the matrix left after
removing all but the first two rows and removing all but the first two columns - here we
have the original matrix H which is a 2 x 2 matrix.

H2 : LPM (H, 2);

H1 and H2 are both positive, so we have a relative minimum. d²f is a positive definite
quadratic form at the given critical point

### plotCP (expr, critPt)

The Maxima function plotCP (f, critPt) is defined in Econ1.mac, and is useful for
the n = 2 case of unconstrained minima or unconstrained maxima. We used this function
in Dowling, Ch. 5 (Dowling05.wxmx). This function shows a 3d plot of the surface f(x,y) 
near the given critical point.

[f, cp1];

plotCP (f, cp1)$

1. Left click the plot once to select it; a border will appear. 
2. Right click once to get a menu.
3. Left click once on "Popout interactively". 
4. Expand the plot to full screen.
5. Rotate the plot using your cursor (holding left button down), to see surface behavior 
    from different angles.
6. Close the gnuplot window to return to the wxMaxima worksheet.

### CPtest (expr, critPts)

The case of the nature of the candidate extremum point (or list of points) of an
expression or function is also evaluated by the function CPtest (expr, critPts), 
defined in Econ1.mac. 

CPtest  internally uses the Hessian matrix method, looking at the leading principal minors.
If there are only two variables in the expression expr, CPtest looks at the signs of 
both zxx and  zyy, and is able to distinguish an inflection point from a saddle point.

CPtest (f, cp1)$

### eigenvalues (numerical-Hessian-matrix)

We can also use the eigenvalues test for a minimum: if all eigenvalues are positive,
we have a relative minimum

We use the Maxima function eigenvalues (numerical-Hessian-matrix), which, if successful, 
returns a list of two lists, the first being a list of eigenvalues found, the second being a list 
of the multiplicities of each of the eigenvalues found.

[eivals, eimult] : eigenvalues (H);

We are interested in the numerical values of the eigenvalues.

float (eivals);

or, you can use (in interactive work)

eivals, numer;

Since all eigenvalues of the numerical Hessian matrix (evaluated at the critical point) are 
positive, the Maxima expression f(x,y) has a relative minimum at the critical point.

### Analyze (expr, critPt)

The function Analyze (expr, critpts) is limited to dealing with an expression which is a 
function of just two variables (as we have here) (see also our Ch. 5 work):

Analyze (f, cp1)$

### Ch. 12, Example 4: Extrema of f (x1, x2, x3) 

Find the critical point(s) and their nature for f(x1,x2,x3) defined by the expression:
     f = -5*x1^2 + 10*x1 + x1*x3 - 2*x2^2 + 4*x2 + 2*x2*x3 - 4*x3^2.

First the quick way with our Maxima function optima (f,vL).

varList : [x1, x2, x3];
f : -5*x1^2 + 10*x1 + x1*x3 - 2*x2^2 + 4*x2 + 2*x2*x3 - 4*x3^2;
optimum (f, varList);

cp;

float(cp);

Next take the longer path, starting with the first derivatives of f.

gradf : jacobian ([f], varList )[1];

J : jacobian ( gradf, varList );

JD : determinant (J);

Since the Jacobian determinant is NOT equal to zero, the three first derivative equations
diff (f,xk) = 0, are functionally independent, and we expect to get a solution from solve.

solns : solve ( gradf, varList);

cp : solns [1];

cp, numer;

Value of f at cp:

float (at (f, cp));

Hessian matrix

H : hessian (f, varList);

The Hessian matrix H is purely numerical in this example, so we skip the normal step of
evaluating H at the critical point cp found, and proceed to the next step of calculating the
leading principal minors of H as it stands.

n = # of variables = 3, so we evaluate the three leading principal minors of the numerical
matrix H.

[LPM (H,1), LPM (H, 2), LPM (H, 3) ];

The sign pattern of the three leading principal minors implies the Taylor series for d²f,
evaluated at the critical point, is a negative definite quadratic in (dx1, dx2, dx3) which 
implies a relative maximum at cp.

CPtest (f, cp);

We can make a 3d plot of two variables at a time near the critical point
using plotCP (expr, critPt). One way to carry this out is to define a Maxima
function F(x1, x2, x3) and call plotCP using this Maxima function with one of the
arguments replaced by one of the three coordinates of the critical point.

cpn : float (cp);

[x10, x20, x30] : map ('rhs, cpn);

We need to place two single quotes (') in front of the expression symbol f to force 
evaluation in this definition of the Maxima function F.

F (x1, x2, x3) := ''f;

Let cpn12 be a two variable critical point list for the variables (x1, x2).

cpn12 : rest (cpn, -1);

In the following call to plotCP, we replace x3 by the number x30 = 0.43...

plotCP (F (x1, x2, x30), cpn12);

Next we hold x2 constant at x20.

cp13 : [cpn[1], cpn[3] ];

plotCP (F (x1,x20, x3), cp13)$

Finally, we hold x1 fixed at x10:

cp23 : rest (cpn, 1);

plotCP (F (x10, x2, x3), cp23)$

### Problem 12.8  Extrema of f(x1, x2, x3) 

Investigate optimizing critical points for
f = 3*x1^2 - 5*x1 - x1*x2 + 6*x2^2 - 4*x2 + 2*x2*x3 + 4*x3^2 + 2*x3 - 3*x1*x3.

varList : [x1, x2, x3];
f : 3*x1^2 - 5*x1 - x1*x2 + 6*x2^2 - 4*x2 + 2*x2*x3 + 4*x3^2 + 2*x3 - 3*x1*x3;
optimum (f, varList);

Dowling's solution (in text, by hand) finds one critical point, for which x3 < 0, which
optimum ignores. Recall that optimum is designed for Economists, who deal with
variables which cannot be negative (usually). A version called optimumAll, with the
same syntax as optimum, does not apply the non-negative filter.

### optimumAll (f, varL) Use Example

optimum is really designed for Economics problems for which the 
search variables are assumed to be non-negative. 
Use optimumAll (f, varList) to accept ALL solutions found by solve.

optimumAll (f, varList);

float (cp);

at (f, %[1]);

### Finding and Evaluating Critical Points "by hand"...

Now for the long way:

gradf : jacobian ([f], varList )[1];

determinant (jacobian (gradf, varList ) );

The Jacobian determinant is not zero, so the three equations implied by the list
gradf, ie., gradf = [0,0,0], are functionally independent, and we expect solve to find 
a symbolic solution.

varList;

solns : solve (gradf, varList);

cp : solns[1];

at (f, cp), numer;

CPtest (f, cp);

CPtest uses the leading principal minors of the Hessian matrix evaluated at the
critical point, but we can check in detail:

H : hessian (f, varList);

Since H is purely numerical, we skip the step of evaluating H at the critical point cp,
and proceed with calculating the leading principal minors of this 3 x 3 matrix.

[LPM (H,1), LPM (H, 2), LPM (H,3)];

All three leading principal minors are positive, so the Taylor expansion d²f about
the critical point is a positive definite quadratic in (dx1, dx2, dx3) and the expression
f has a relative minimum at that point and f is locally convex.

As in Example 4 above, we can call plotCP for 3d views of this minimum after some prep:

[x10, x20, x30] : float (map ('rhs, cp));

Use two single quotes (') in front of the expression f to define the Maxima function F.

F( x1, x2, x3) := ''f;

grind(%)$

For fixed x3 = x30, look at (x1,x2) local behavior:

plotCP (F (x1, x2, x30), rest (cp, -1))$

For fixed x2 = x20, look at (x1, x3) behavior:

plotCP (F (x1, x20, x3), [cp[1], cp[3]] )$

For fixed x1 = x10, look at (x2, x3) behvior:

plotCP (F (x10, x2, x3), [cp[2], cp[3]] )$