# Which optimization algorithm converges faster?

everyone. I'm having a large scale unconstrained optimization problem. If I treat the unconstrained problem as a constrained problem with infinity constraints, I should be able to use both the fminunc and fmincon function in MATLAB. When using the fminunc function, I should provide the gradient and the sparse pattern of the Hessian. While using the fmincon function, I can choose the l-bfgs method to approximate the Hessian.

So, here is my question. Which one converges faster? Are there any good reasons to choose one over another? Or any principles to guide such choice?

Update: Actually the problem is a least-squares problem. However the number of terms is big, so the Jacobian would be too large to store. This is the reason why I don't use the lsqnonlin function. The algorithm used in fminunc for large scale problem is a trust-region method (details can be found in fminunc documentation), and the algorithm in fmincon is l-bfgs (see fmincon documentation). What I'm asking is more of a generative comparison because there are many C/C++ implementations of these algorithms. Answers may not be necessarily constrained by MATLAB. Thanks!

Update 2: I was unclear, the form of my problem is like this: $$\min_W (\sum_i\|Y_i-WX_i\|_2^2).$$

Update 3: Thank you all! littleO is right. But I'm still curious about the efficiency of the methods mentioned above under general condition.

• My suggestion is to run a minimal test problem that inherits the structure of a larger problem you wanna solve. Otherwise it is hard to tell, I would put my bet on fminunc. Aug 6, 2013 at 3:37
• Since your problem is unconstrained, I'd guess fminunc would be faster. But can you tell us more about the function you're minimizing? Aug 6, 2013 at 4:48
• Are $Y_i$ and $X_i$ matrices? If so, are you using $\| \cdot \|_2$ to denote what is usually called the Frobenius norm of a matrix? Or is it actually an induced norm ? Aug 6, 2013 at 8:04
• I think this minimization problem has a closed-form solution! Aug 6, 2013 at 8:06
• Since you're using a Frobenius norm, it will be possible to express your problem as a least squares problem, and then solve it using a least squares solver from numerical linear algebra. This should be much more efficient than using a general purpose solver. Aug 6, 2013 at 8:10

If your problem is truly a least squares problem: $$\text{minimize}_x \quad \|Ax - b\|_2^2$$ then you should probably use a technique from numerical linear algebra, rather than some general purpose optimization algorithm. Least squares problems have been studied extensively in numerical linear algebra and very good algorithms / libraries exist for least squares problems.

In Matlab, you can solve least squares problems just using the backslash operator:

x = A\b ;

$A$ can be stored as a sparse matrix in Matlab if it's huge.

Edit: Your specific problem can be expressed as a least squares problem. First note that \begin{align} \|Y_i - W X_i\|_F &= \| Y_i^T - X_i^T W^T \|_F \\ &= \| B^i - A^i Z \|_F \end{align} where I have changed notation in the last line. (Note that $W^T = Z$.)

Furthermore, \begin{align} \| B^i - A^i Z \|_F &= \left \| \begin{bmatrix} B^i_1 - A^i Z_1 \\ \vdots \\ B^i_N - A^i Z_N \end{bmatrix} \right \|_2 \end{align} where $B^i_1,\ldots,B^i_N$ are the columns of $B^i$ and $Z_1,\ldots,Z_N$ are the columns of $Z$. Let $$\tilde{B}^i = \begin{bmatrix} B^i_1 \\ \vdots \\ B^i_N \end{bmatrix}$$ and $$\tilde{A}^i = \begin{bmatrix} A^i & & \\ & \ddots & \\ & & A^i \end{bmatrix}$$ and $$\tilde{Z} = \begin{bmatrix} Z_1 \\ \vdots \\ Z_N \end{bmatrix}.$$ Then $$\|B^i - A^i Z \|_F = \|\tilde{A}^i \tilde{Z} - \tilde{B}^i \|_2.$$ Note that $\tilde{Z}$ and $\tilde{B}^i$ are column vectors while $\tilde{A}^i$ is a matrix, so we are close to expressing our problem as a standard least squares problem.

Next note that \begin{align} \sum_{i=1}^K \|Y_i - W X_i \|_F^2 &= \sum_{i=1}^K \|\tilde{A}^i \tilde{Z} - \tilde{B}^i \|_2^2. \end{align} Let $$A = \begin{bmatrix} \tilde{A}^1 & & \\ & \ddots & \\ & & \tilde{A}^K \end{bmatrix}$$ and $$b = \begin{bmatrix} \tilde{B}^1 \\ \vdots \\ \tilde{B}^K \end{bmatrix}.$$ We finally see that your problem is equivalent to $$\text{minimize}_{\tilde{Z}} \quad \|A \tilde{Z} - b \|_2^2.$$

An analytical solution is given by $$\tilde{Z} = (A^T A)^{-1} A^T b.$$ ($\tilde{Z}$ satisfies the normal equations.) However, in practice you should probably compute $\tilde{Z}$ using a least squares solver. In Matlab, it's just

Ztilde = A\b ;

• I was unclear. I've added details about my problem. Aug 6, 2013 at 8:02

If you can approximate the action of the Jacobian, I'd try fminunc with the HessMult option. This amounts to trust-region Gauss-Newton. In terms of data storage, you don't need the actual Jacobian of the map (or Hessian of the objective function), only a function that can apply it to vectors. L-BFGS is basically a clever way to approximate the action of the inverse Hessian, so why not actually use the Hessian itself?

In detail, suppose you want to solve the following least-squares optimization problem, $$\min_x \frac{1}{2}||F(x)-y||^2,$$

where $F$ has a Taylor expansion at $p$, $$F(x) = F_0 + J (x-p) + \dots .$$

Then the objective function can be expanded to, $$\frac{1}{2}||F(x)-y||^2 = \frac{1}{2}[F_0 + J (x-p) + \dots -y]^T[F_0 + J (x-p) + \dots -y] \\ =\frac{1}{2}||F_0-y||^2 + (F_0-y)^T J (x-p) + \frac{1}{2}(x-p)^T J^T J(x-p) + \dots \\ = A + G (x-p) + (x-p)^T H (x-p) + \dots .$$

This is not a true second order Taylor series of the objective function because of possible interactions of $A,G,H$ with terms in the $\dots$'s. However, it is the Gauss-Newton approximation which is generally well accepted. The motivation for this is to note that at the true solution, we expect the misfit will be small. Ie, $(F_0-y)^T = O(x-p)"$.

You need functions, say myGradient(p) and myHessian(p,w), that compute $G = (F_0-y)^T J$ and the action $Hw = J^TJw$ on an input vector $w$. Then, simply pass the function handles @myGradient & @myHessian as arguments to fminunc.

fminunc will then apply Newton's method to the optimization problem, $$x_{k+1}=x_k - H^{-1} G x_k.$$

The inverse $H^{-1} G x_k$ is applied through conjugate gradient, each step of which involves evaluating products like $Hw$. Whenever the action of the Hessian, $Hw$, is needed, fminunc will call your myHessian(p,w) function.

Thus you never need to store the full Jacobian of the forward map (Hessian of the objective function).

Edit: This method is exact after 1 Newton step for the linear least squares problem in the updated version of the question, as it amounts to using conjugate gradient on the normal equations.