# Iterative solution for a minimization problem involving matrix equations

I have a real valued function $$F$$ for which I am looking to find its global minimum. The function is well behaved and I can obtain its Jacobian. I could also compute the Hessian but the function depends on a large number of variables so it's not practical. We can still assume that the Hessian is definite positive over the whole domain so there is only one minimum (provided that there is a point where the Jacobian is null).

We can write the function $$F$$ to depend on two matrices, $$A$$ and $$B$$, where $$A$$ is a definite positive matrix and $$B$$ is a matrix with no particular restriction.

I have an analytic form for $$F_A=\frac{\partial}{\partial A}F$$ and $$F_B=\frac{\partial}{\partial B}F$$. The solution I am thus looking for is the values $$A$$ and $$B$$ that solve $$F_A=0$$ and $$F_B=0$$ simultaneously.

I can find the zero for $$F_A(A,B_0)$$ where $$B_0$$ is fixed (it's this equation) and conversely I can find the zero for $$F_B(A_0,B)$$ where $$A_0$$ is fixed (it's a Sylvester equation) but finding both at the same time doesn't seem to be possible analytically. So I am thinking to iterate in the following way:

1. Provide an initial guess $$A_0$$ and $$B_0$$
2. Find the new $$A$$ so $$F_A(A,B_0)=0$$, this becomes the new $$A_0$$
3. Find the new $$B$$ so $$F_B(A_0,B)=0$$, this becomes the new $$B_0$$
4. Iterate until we converge to some solution within some tolerance

Is there a more efficient approach to root finding for such problem? Computing the Hessian - so I could use a Newton-Raphson approach - might be doable but the Hessian dimension will be massive.

• Perhaps consider using a gradient based quasi-Newton optimization approach such as BFGS to find the minimum? At the minimum you should have $F_A=0$ and $F_B=0$. Since you say you can assume positive definiteness of the Hessian, BFGS should get to the global minimum.
– NNN
Jul 9, 2022 at 4:28
• If your function F depends on two matrices A(n,n) and B(m,m), doesn't this just mean that it is a function of $n^2$ parameters $A_{i,j}$ and $m^2$ parameters $B_{i,j}$? So we need to search for a minimum of F in the space of those parameters. Certainly there are methods of minimization that don't require the Hessian (and some don't even need the Jacobian), but if the space is highly dimensional then this may become a difficult problem (depending on the function F). I would look in the literature on protein folding, they solve it as a global minimization problem with hundreds of variables. Jul 9, 2022 at 4:50
• Methods like the one suggested tend to converge slowly, but convergence can be accelerated significantly using methods such as Anderson acceleration. See, for example, eprints.ma.man.ac.uk/2360 Jul 11, 2022 at 4:13
• @AmitHochman yes, I do notice that the solution seems to "zigzag" towards the minimum, and this is not fast.
– PC1
Jul 11, 2022 at 5:27

$$\def\a{\alpha}\def\b{\beta}\def\g{\gamma}\def\t{\theta} \def\l{\lambda}\def\s{\sigma}\def\e{\varepsilon} \def\n{\nabla}\def\o{{\tt1}}\def\p{\partial} \def\E{{\cal E}}\def\F{{\cal F}}\def\G{{\cal G}} \def\LR#1{\left(#1\right)} \def\BR#1{\Big(#1\Big)} \def\bR#1{\big(#1\big)} \def\qiq{\quad\implies\quad} \def\vc#1{\operatorname{vec}\LR{#1}} \def\trace#1{\operatorname{Tr}\LR{#1}} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\c#1{\color{red}{#1}} \def\m#1{\left[\begin{array}{r}#1\end{array}\right]} \def\dk#1{\LR{#1_k-#1_{k-\o}}}$$Let's use uppercase to denote matrix quantities, lowercase to denote vectors, and Greek letters to denote scalar quantities. Further, we'll reserve superscripts for conjugations (or transpositions) and subscripts for components. This convention requires renaming some of the problem variables \eqalign{ \bR{F,F_A,F_B} \to \bR{\phi,P,Q} \\ } The matrix inner product is denoted by a colon \eqalign{ A:P &= \sum_{i=1}^m\sum_{j=1}^n A_{ij}P_{ij} \;=\; \trace{A^TP} \\ A:A &= \|A\|^2_F \\ } Another convenient notation, is that vectorizing a matrix creates a vector of the same name \eqalign{ a = \vc A,\quad b = \vc B,\quad etc. \\ } One last point is that the constrained matrix $$A$$ needs to be eliminated in favor of an unconstrained matrix variable $$C,\,$$ i.e.
\eqalign{ A &= CC^T \qiq dA = \LR{dC\,C^T +\, C\,dC^T} \\ P:dA &= P:\LR{dC\,C^T +\, C\,dC^T} \\ &= \LR{P+P^T}:dC\,C^T \\ &= 2PC:dC \\ &\doteq R:dC \\ }
Since the gradients of the function are known, the differential can be written and vectorized \eqalign{ d\phi &= P:dA \;+\; Q:dB \\ &= R:dC \;+\; Q:dB \\ &= \vc{R}:\vc{dC} \;+\; \vc{Q}:\vc{dB} \\ &= r:dc \;+\; q:db \\ &= \m{q\\r}:\m{db\\dc} \\ &= g:dx \\ } where the last line introduces the partitioned vectors \eqalign{ x &= \m{b \\ c}=\m{\vc{B} \\ \vc{C}}, \qquad g = \m{q \\ r} = \grad{\phi}{x} \\ } The reason for writing the problem in this form is that you now have an optimization problem in terms of a single unconstrained vector variable, so you can employ any standard gradient-based method for its solution.
Initialize \eqalign{ x_0 &= random \qquad\qquad\qquad\qquad\qquad\quad \\ } First step \eqalign{ g_0 &= g(x_0) \\ \phi_0 &= \phi(x_0) \\ x_1 &= x_0 - \LR{\frac{0.05\;\big|\phi_0\big|}{g_0^Tg_0}} g_0 \qquad\qquad\quad \\ k &= \o \\ } Subsequent steps \eqalign{ g_k &= g(x_k) \\ x_{k+1} &= x_k - \LR{\frac{\dk{x}^T\dk{g}}{\dk{g}^T\dk{g}}}g_k \\ k &= k+\o \\ } Stop when $$g_k\approx 0$$