[Math] Matlab: Gradient and Hessian of a function with vector input of user specified size

hessian-matrixMATLABnonlinear optimizationoptimizationvector analysis

I need to write a matlab m file that takes the following function of $x=(x_{1},x_{2},\cdots, x_{2n})$, and for $n=10$, $n=100$, $n=500$, $$f(x) = \frac{1}{2}\sum_{i=1}^{2n}i(x_{i})^{2}-\sum_{i=1}^{2n}x_{i}+\sum_{i=2}^{2n}\left[\frac{1}{4}\left(x_{i}+x_{i-1} \right)^{2} +\left(x_{i}-x_{i-1}\right)^{2}\right]+\left(x_{2n}-x_{1} \right)^{2} $$ calculates the gradient, and each component $H_{ij}=\frac{\partial^{2}f}{\partial x_{i} \partial x_{j}}$ of the Hessian matrix.

Now, the course that I am doing this for is NOT a course on Matlab (its a nonlinear optimization course). In fact, they take for granted that we all are extremely proficient in Matlab and are comfortable programming to this sophisticated of a degree, which I am not.

Therefore, I am asking, based on the code that I have done thus far, which I am including in a grapic, how do I write the code to

Calculate the gradients, based on the number inputted into the function by the user, and put all the components into a function of the form $g = [\cdots; \cdots; \cdots]$ with $m$ components.
Calculate the Hessian components, also based on the number inputted into the function by the user. EDIT: I think I might have come up with suitable code for this. Please see my screenshot and tell me if it looks okay. (Still have no idea how to do this for the gradient though)

I am assuming I would need to set up some kind of dynamic vector for each of these things to be inputted into, and I have absolutely no idea how to do this. So, any assistance with which you could provide me would be greatly beneficial.

Thank you for your time and patience.

Best Answer

Your current solution is based on symbolic math. There are two ways to move forward.

The first is to have continue with symbolic math and use automated differentiation. That will not be the fastest solution, because for $n=500$ you have roughly half a million elements in the Hessian you need to compute.

The second is to observe that your function is quadratic in $x$, and can be written as $0.5 x^TAx + b^Tx + c$. After computing $A$ and $b$, the gradient can be computed efficiently as $Ax+b$, and the Hessian is $A$.

Your objective function is: $$f(x) = \frac{1}{2}\sum_{i=1}^{2n}i(x_{i})^{2}-\sum_{i=1}^{2n}x_{i}+\sum_{i=2}^{2n}\left[\frac{1}{4}\left(x_{i}+x_{i-1} \right)^{2} +\left(x_{i}-x_{i-1}\right)^{2}\right]+\left(x_{2n}-x_{1} \right)^{2}. $$ The easiest way to set up your code is to initialize $A$ and $b$ to a matrix or vector with all zeros:

A = zeros(2*n,2*n);
b = zeros(2*n, 1);

and to go through the terms of $f$ and update $A$ and $b$. For example, the first term affects the diagonal of $A$:

for i = 1:2*n
  A(i,i) = A(i,i) + i;
end

You should be able to figure out the code for the other terms. Note that you should compute $A$ and $b$ only once, not repeatedly in each iteration.

Related Solutions

[Math] Write a function in MATLAB that returns the value of a function f (where f needs to be written by the user)

So you want to write a subroutine that returns the jacobian matrix at point $x\in \mathbb{R}^2$ of a function handle @f in matlab. Also you wanna do this numerically, not symbolically.

Here is a simple finite difference subroutine that does the job for f being a function handle:

function J=numjacobian(f,x)

n=length(x); % n=2 in your case, but can be higher 
fx=feval(f,x); % evaluate the function at point x
step=1e-6; % difference step you choose, can be 1e-10 if you like

for i=1:n
   xstep = x;
   % compute the i-th component partial derivative 
   % numerically using first order forward difference approximation
   xstep(i)=x(i)+step;
   J(:,i)=(feval(f,xstep)-fx)/step;
end;
end

The input f is a function handle, x is any point. The output J is the jacobian matrix at x. Example: $f = (x^2+y^2, x^2-y^2)$:

  >>f = @(x) [x(:,1).^2 + x(:,2).^2, x(:,1).^2 - x(:,2).^2]; 
  >>numjacobian(f,[1 2])
    ans =
      2.00000100036846          4.00000100064801
      2.00000099992437         -4.00000100064801

For the Symbolic Toolbox approach, please refer to my answer here in this question: Solving a system with Newton's method in matlab? . A more object-oriented approach would be like you mentioned, allowing user to input any type of function: inline function, symbolic, function handle, etc. Converting symbolic function to function handle can be done by matlabFunction in MATLAB. And inline functions are string, hence easy to be converted to a symbolic function.

[Math] How to find gradient and hessian of summation

Once you have convinced yourself that $$\frac{\partial f(x)}{\partial x_i}=\frac{\partial \|x\|}{\partial x_i}=\frac{x_i}{\|x\|}$$ Then, recall that, the Hessian matrix of $f$ is the $n\times n$ matrix $\textbf H$ with the $(i,j)$-entry given by $$\textbf{H}_{ij}=\frac{\partial^2 f(x)}{\partial x_j\partial x_i}=\frac{\partial}{\partial x_j}\bigg(\frac{\partial f(x)}{\partial x_i}\bigg)$$ thus, we will have to calculate the latter in order to give the general input of the matrix.

Using the quotient rule, we see that $$\begin{align} \frac{\partial}{\partial x_j}\bigg(\frac{\partial f(x)}{\partial x_i}\bigg) &= \frac{\partial}{\partial x_j}\bigg(\frac{x_i}{\|x\|}\bigg)=\cfrac{\cfrac{\partial x_i}{\partial x_j}\cdot\|x\|-x_i\cdot\cfrac{\partial \|x\|}{\partial x_j}}{\|x\|^2} \\ &= \frac{1}{\|x\|}\frac{\partial x_i}{\partial x_j}-\frac{x_ix_j}{\|x\|^3} \end{align}$$ Note that, $$\frac{\partial x_i}{\partial x_j} = \delta_{ij} =\begin{cases} 1 & \textrm{if } i=j \\ 0 & \textrm{if } i\neq j \end{cases}$$ Hence, we have $$\mathbf{H}_{ij}=\frac{1}{\|x\|}\delta_{ij}-\frac{x_ix_j}{\|x\|^3} = \frac{\delta_{ij}\|x\|^2-x_ix_j}{\|x\|^3}$$ as, for example, if $n=2$, then $f$ is given by $f(x,y)=\|(x,y)\|=\sqrt{x^2+y^2}$ and then $$\mathbf{H}=\begin{pmatrix} \cfrac{\delta_{11}\|(x,y)\|^2-x^2}{\|(x,y)\|^3} & \cfrac{\delta_{12}\|(x,y)\|^2-xy}{\|(x,y)\|^3} \\ \cfrac{\delta_{21}\|(x,y)\|^2-yx}{\|(x,y)\|^3} & \cfrac{\delta_{22}\|(x,y)\|^2-y^2}{\|(x,y)\|^3} \end{pmatrix}=\begin{pmatrix} \cfrac{y^2}{(x^2+y^2)^{3/2}} & -\cfrac{xy}{(x^2+y^2)^{3/2}} \\ -\cfrac{xy}{(x^2+y^2)^{3/2}} & \cfrac{x^2}{(x^2+y^2)^{3/2}} \end{pmatrix} = \frac{1}{(x^2+y^2)^{3/2}}\begin{pmatrix} y^2&-xy \\ -xy&x^2 \end{pmatrix}$$

Best Answer

Related Solutions

[Math] Write a function in MATLAB that returns the value of a function f (where f needs to be written by the user)

[Math] How to find gradient and hessian of summation

Related Question