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

Question

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

1. 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.
2. 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.

Answer 1

Hint.

Taking a symbolic processor we can build and generalize the tables for $\nabla f$ and $\nabla^2 f$

With MATHEMATICA we have for $n = 6$

n = 6;
X = Table[Subscript[x, k], {k, 1, 2 n}];
f = 1/2 Sum[k Subscript[x, k]^2, {k, 1, 2 n}] - Sum[Subscript[x, k], {k, 1, 2 n}] + Sum[1/4 (Subscript[x, k] + Subscript[x, k - 1])^2 + (Subscript[x, k] - Subscript[x, k - 1])^2, {k, 2, 2 n}] + (Subscript[x, 2 n] - Subscript[x, 1])^2;
grad = Grad[f, X];
H = Grad[grad, X];
MatrixForm[2 grad]
MatrixForm[2 H]

so we have

$$ \nabla f = \frac 12 \left( \begin{array}{c} -3 x_2+11 x_1-4 x_{12}-2 \\ -3 x_1+14 x_2-3 x_3-2 \\ -3 x_2+16 x_3-3 x_4-2 \\ -3 x_3+18 x_4-3 x_5-2 \\ -3 x_4+20 x_5-3 x_6-2 \\ -3 x_5+22 x_6-3 x_7-2 \\ -3 x_6+24 x_7-3 x_8-2 \\ -3 x_7+26 x_8-3 x_9-2 \\ -3 x_8+28 x_9-3 x_{10}-2 \\ -3 x_9+30 x_{10}-3 x_{11}-2 \\ -3 x_{10}+32 x_{11}-3 x_{12}-2 \\ -3 x_{11}+33 x_{12}-4 x_1-2 \\ \end{array} \right) $$

$$ \nabla ^2 f = \frac 12\left( \begin{array}{cccccccccccc} 11 & -3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -4 \\ -3 & 14 & -3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -3 & 16 & -3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -3 & 18 & -3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & -3 & 20 & -3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & -3 & 22 & -3 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & -3 & 24 & -3 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & -3 & 26 & -3 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & -3 & 28 & -3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -3 & 30 & -3 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -3 & 32 & -3 \\ -4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -3 & 33 \\ \end{array} \right) $$

Those tables can easily be generalized for any $n$ so in the programmed function those generalized tables can be directly inserted. As can be observed $H$ is a sparse diagonal dominant matrix that can be conveniently compacted before handling.

And for any $n$

$$ \nabla f = \frac 12 \left( \begin{array}{c} -3 x_2+11 x_1-4 x_{2n}-2 \\ -3 x_1+14 x_2-3 x_3-2 \\ -3 x_2+16 x_3-3 x_4-2 \\ -3 x_3+18 x_4-3 x_5-2 \\ -3 x_4+20 x_5-3 x_6-2 \\ \vdots\\ -3 x_{2n-4}+(2(2n-1)+11-4) x_{2n-3}-3 x_{2n-2}-2 \\ -3 x_{2n-3}+(2(2n-1)+11-2) x_{10}-3 x_{2n-1}-2 \\ -3 x_{2n-2}+(2(2n-1)+11-1) x_{2n-1}-3 x_{2n}-2 \\ -3 x_{2n-1}+(2(2n-1)+11) x_{2n}-4 x_1-2 \\ \end{array} \right) $$

$$ \nabla ^2 f = \frac 12\left( \begin{array}{cccccccccccc} 11 & -3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -4 \\ -3 & 14 & -3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & -3 & 16 & -3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & -3 & 18 & -3 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots\\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -3 & 4n+6 & -3 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -3 & 4n+8 & -3 \\ -4 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & -3 & 2(2n-1)+11 \\ \end{array} \right) $$

NOTE

The matlab function follows

function [val,g,H] = nTermSumFunction(x,n)
  g  = zeros(2*n,1);
  H  = zeros(2*n,2*n);
  k0 = 14;

  g(1)     =(-3*x(2)+11*x(1)-4*x(2*n)-2)/2;
  g(2*n)   =(-3*x(2*n-1)+(4*n+9)*x(2*n)-2)/2;
  H(1,1)   = 11/2; H(1,2) = -3/2; H(1,2*n) = -4/2;
  H(2*n,1) = -4/2; H(2*n,2*n-1) = -2/2;H(2*n,2*n) = (4*n+9)/2;
  val      = (x(2*n)-x(1))^2+1/2*x(1)^2-x(1);

  for k = 2:2*n-1
    val      = val + k/2*x(k)^2-x(k)+(x(k)+x(k-1))^2/4+(x(k)-x(k-1))^2;
    g(k)     = (-3*x(k-1)+k0*x(k)-3*x(k+1)-2)/2;
    H(k,k)   = k0/2;
    H(k,k+1) = -3/2;
    H(k,k-1) = -3/2;
    k0       = k0 + 2;
  end

Answer 2

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.