Main Content

Find Extremum of Multivariate Function and Its Approximation

This example shows how to find the extremum of a multivariate function and its approximation near the extremum point. This example uses symbolic matrix variables to represent the multivariate function and its derivatives. Symbolic matrix variables are available starting in R2021a.

Consider the multivariate function f(x)=sin(xTAx), where x is a 2-by-1 vector and A is a 2-by-2 matrix. To find a local extremum of this function, compute a root of the derivative of f(x). In other words, find the solution of the derivative f(x0)=0.

Create Function and Find Its Derivative

Create the vector x and the matrix A as symbolic matrix variables. Define the function f(x)=sin(xTAx).

syms x [2 1] matrix
syms A [2 2] matrix
f = sin(x.'*A*x)
f = sin(xTAx)sin(transpose(symmatrix('x', [2 1]))*symmatrix('A', [2 2])*symmatrix('x', [2 1]))

Compute the derivative D of the function f(x) with respect to the vector x. The derivative D is displayed in compact matrix notation in terms of x and A.

D = diff(f,x)
D = cos(xTAx)xTA+xTATcos(transpose(symmatrix('x', [2 1]))*symmatrix('A', [2 2])*symmatrix('x', [2 1]))*(transpose(symmatrix('x', [2 1]))*symmatrix('A', [2 2]) + transpose(symmatrix('x', [2 1]))*transpose(symmatrix('A', [2 2])))

Convert symmatrix Objects to sym Objects

The symbolic matrix variables x, A, f, and D are symmatrix objects. These objects represent matrices, vectors, and scalars in compact matrix notation. To show the components of these variables, convert the symmatrix objects to sym objects using symmatrix2sym.

xsym = symmatrix2sym(x)
xsym = 

(x1x2)[x1; x2]

Asym = symmatrix2sym(A)
Asym = 

(A1,1A1,2A2,1A2,2)[A1_1, A1_2; A2_1, A2_2]

fsym = symmatrix2sym(f)
fsym = (sin(x1A1,1x1+A1,2x2+x2A2,1x1+A2,2x2))[sin(x1*(A1_1*x1 + A1_2*x2) + x2*(A2_1*x1 + A2_2*x2))]
Dsym = symmatrix2sym(D)
Dsym = (cos(x1A1,1x1+A1,2x2+x2A2,1x1+A2,2x2)2A1,1x1+A1,2x2+A2,1x2cos(x1A1,1x1+A1,2x2+x2A2,1x1+A2,2x2)A1,2x1+A2,1x1+2A2,2x2)[cos(x1*(A1_1*x1 + A1_2*x2) + x2*(A2_1*x1 + A2_2*x2))*(2*A1_1*x1 + A1_2*x2 + A2_1*x2), cos(x1*(A1_1*x1 + A1_2*x2) + x2*(A2_1*x1 + A2_2*x2))*(A1_2*x1 + A2_1*x1 + 2*A2_2*x2)]

Substitute Numeric Values and Find the Minimum

Suppose you are interested in the case where the value of A is [2 -1; 0 3]. Substitute this value into the function fsym.

fsym = subs(fsym,Asym,[2 -1; 0 3])
fsym = sin(3x22+x12x1-x2)sin(3*x2^2 + x1*(2*x1 - x2))

Substitute the value of A into the derivative Dsym

Dsym = subs(Dsym,Asym,[2 -1; 0 3])
Dsym = (cos(3x22+x12x1-x2)4x1-x2-cos(3x22+x12x1-x2)x1-6x2)[cos(3*x2^2 + x1*(2*x1 - x2))*(4*x1 - x2), -cos(3*x2^2 + x1*(2*x1 - x2))*(x1 - 6*x2)]

Then, apply the symbolic function solve to get a root of the derivative.

[xmin,ymin] = solve(Dsym,xsym,'PrincipalValue',true);
x0 = [xmin; ymin]
x0 = 

(00)[sym(0); sym(0)]

Plot the function f(x) together with the extremum solution x0. Set the plot interval to -1<x1<1 and -1<x2<1 as the second argument of fsurf. Use fplot3 to plot the coordinates of the extremum solution.

fsurf(fsym,[-1 1 -1 1])
hold on
fplot3(xmin,ymin,subs(fsym,xsym,x0),'ro')
view([-68 13])

Figure contains an axes. The axes contains 2 objects of type functionsurface, parameterizedfunctionline.

Approximate Function Near Its Minimum

You can approximate a multivariate function around a point x0 with a multinomial using the Taylor expansion.

f(x)f(x0)+f(x0)(x-x0)+12(x-x0)TH(f(x0))(x-x0)

Here, the term f(x0) is the gradient vector, and H(f(x0)) is the Hessian matrix of the multivariate function f(x) calculated at x0.

Find the Hessian matrix and return the result as a symbolic matrix variable.

H = diff(f,x,x.')
H = -sin(xTAx)ATx+AxxTA+xTAT+cos(xTAx)AT+A- sin(transpose(symmatrix('x', [2 1]))*symmatrix('A', [2 2])*symmatrix('x', [2 1]))*(transpose(symmatrix('A', [2 2]))*symmatrix('x', [2 1]) + symmatrix('A', [2 2])*symmatrix('x', [2 1]))*(transpose(symmatrix('x', [2 1]))*symmatrix('A', [2 2]) + transpose(symmatrix('x', [2 1]))*transpose(symmatrix('A', [2 2]))) + cos(transpose(symmatrix('x', [2 1]))*symmatrix('A', [2 2])*symmatrix('x', [2 1]))*(transpose(symmatrix('A', [2 2])) + symmatrix('A', [2 2]))

Convert the Hessian matrix H(f(x0)) to the sym data type, which represents the matrix in its component form. Use subs to evaluate the Hessian matrix for A = [2 -1; 0 3] at the minimum point x0.

Hsym = symmatrix2sym(H);
Hsym = subs(Hsym,Asym,[2 -1; 0 3]);
H0 = subs(Hsym,xsym,x0)
H0 = 

(4-1-16)[sym(4), -sym(1); -sym(1), sym(6)]

Evaluate the gradient vector f(x0) at x0.

D0 = subs(Dsym,xsym,x0)
D0 = (00)[sym(0), sym(0)]

Compute the Taylor approximation to the function near its minimum.

fapprox = subs(fsym,xsym,x0) + D0*(xsym-x0) + 1/2*(xsym-x0).'*H0*(xsym-x0)
fapprox = 

x12x1-x22-x2x12-3x2x1*(2*x1 - x2/2) - x2*(x1/2 - 3*x2)

Plot the function approximation on the same graph that shows f(x) and x0.

hold on
fsurf(fapprox,[-1 1 -1 1])
zlim([-1 3])

Figure contains an axes. The axes contains 3 objects of type functionsurface, parameterizedfunctionline.