%% Unconstrained minimization of quadratic forms
% The quadratic function to be minimized is parameterized by a matrix Q and a vector c
% f(\mathbf x) = \frac{1}{2}\mathbf x^T\mathbf Q\mathbf x + \mathbf c^T\mathbf x
% Below we analyze the function for a few particular matrices Q (and vectors r).
% We find the critical points by solving for x where the gradient vanishes
% \nabla f(\mathbf x) = \mathbf 0
% and examine if the function achieve a minimum or maximum or if the stationary point is just a saddle point by analyzing the spectrum (eigenvalues) of the Hessian (matrix of the second derivatives, which in this case is just Q).
%% Minimum at the critical point
% First, the data for the problem

Q = [1 1; 1 2]; c = [0; 1];

% Now we compute the critical point just by solving the linear system

xst = -Q\c

% Checking if the critical point is indeed a minimum can be done by computing the eigenvalues

eig(Q)

% But for larger problems, computing the eigenvalues is an overkill. Instead, it can be just checked if Cholesky factorization exists

[~,p] = chol(Q)

%
%% Plotting the function

x1 = -10:0.1:10;
x2 = x1;

[X1,X2] = meshgrid(x1,x2);

Z = 1/2*(Q(1,1)*X1.^2 + 2*Q(1,2)*X1.*X2 + Q(2,2)*X2.^2) + c(1)*X1 + c(2)*X2;

[Ca,ha] = contour(X1,X2,Z);
set(gca,'FontSize',12)
set(ha,'ShowText','on','TextStep',get(ha,'LevelStep'))
hold on, plot(xst(1),xst(2),'x'), hold off
xlabel('x_1'), ylabel('x_2')
grid on
surfc(X1,X2,Z)
set(gca,'FontSize',12)
xlabel('x_1'), ylabel('x_2'), zlabel('f(x)')

%
%% Finding the minimizer using Optimization Toolbox –⁠ solver-based workflow

clear x

% See https://www.mathworks.com/help/optim/, https://www.mathworks.com/help/optim/optimization-problem-setup-solver-based.html

f = @(x) 1/2*x'*Q*x + c'*x
x0 = [1,1]';

[xopt,fval] = fminunc(f,x0)

%
%% Finding the minimizer using Optimization toolbox –⁠ problem-based workflow

clear x x0

% See https://www.mathworks.com/help/optim/problem-based-approach.html. The optimization problem is formulated by creating optimproblem object first and then populating it with the objetive function and constraints.

prob = optimproblem('ObjectiveSense','minimize')
x = optimvar('x',2,1,'Type','continuous')
prob.Objective = 1/2*x'*Q*x + c'*x
show(prob)

% Once the problem is defined, it can be solved

x0.x = [1;1]
sol = solve(prob,x0)

%
%% Finding the minimizer using CVX

clear x

% See http://cvxr.com/cvx/

cvx_begin
    variable x(2)
    minimize( 1/2*quad_form(x,Q)+c'*x )
cvx_end

%
%% When critical point turns out a saddle point

Q = [-1 1; 1 2];  c = [0; 1];

xopt = -Q\c
eig(Q)

% Obviously the signs of the two eigenvalues differ, the stationary point is a saddle point.
% That the point is not a minimum can also be checked in a cheaper way by attempting to do Cholesky decomposition of the matrix

[~,p] = chol(Q)

% It fails, which confirms that the matrix Q is not positive definite and the point is not a minimizer.
%
%% Plotting the function

x1 = -10:0.1:10;
x2 = x1;

[X1,X2] = meshgrid(x1,x2);

Z = 1/2*(Q(1,1)*X1.^2 + 2*Q(1,2)*X1.*X2 + Q(2,2)*X2.^2) + c(1)*X1 + c(2)*X2;

[Ca,ha] = contour(X1,X2,Z);
set(gca,'FontSize',12)
set(ha,'ShowText','on','TextStep',get(ha,'LevelStep'))
hold on, plot(xst(1),xst(2),'x'), hold off
xlabel('x_1'), ylabel('x_2')
grid on

surfc(X1,X2,Z)
set(gca,'FontSize',12)
xlabel('x_1'), ylabel('x_2'), zlabel('f(x)')

%
%% Finding the minimizer using Optimization Toolbox
% See https://www.mathworks.com/help/optim/

clear x

f = @(x) 1/2*x'*Q*x + c'*x
x0 = [1,1]';

[xopt,fval] = fminunc(f,x0)

%
%% What if the Hessian is singular?

Q = [2 0; 0 0];
c = [3 0]';

eig(Q)

xst = -Q\c

%
%% Plotting the function

x1 = -10:0.1:10;
x2 = x1;

[X1,X2] = meshgrid(x1,x2);

Z = 1/2*(Q(1,1)*X1.^2 + 2*Q(1,2)*X1.*X2 + Q(2,2)*X2.^2) + c(1)*X1 + c(2)*X2;

[Ca,ha] = contour(X1,X2,Z);
set(gca,'FontSize',12)
set(ha,'ShowText','on','TextStep',get(ha,'LevelStep'))
hold on, plot(xst(1),xst(2),'x'), hold off
xlabel('x_1'), ylabel('x_2')
grid on

surfc(X1,X2,Z)
set(gca,'FontSize',12)
xlabel('x_1'), ylabel('x_2'), zlabel('f(x)')