I need to find the vector that minimizes this matrix equation;

$$\bar{v}.M_1.v-(\bar{v}.M_2.v)^2 $$

$v$ is a normalized complex vector and $\bar{v}$ is the complex conjugate but I can arrange parameters in such a way that it can be all real. $M_1$ and $M_2$ are Hermitian matrices and therefore the function I am asking is eventually a real number.

My idea was to rewrite the vector in terms of coefficients and then calculate the product and I get a quartic polynomial but solving it in higher dimensions is a mess.

So I thought maybe it is possible to use some matrix techniques which always simplifies things. Approximations are also accepted.

Also is there calculus techniques such that I can take derivative of this function and find the optimum vector?

  • $\begingroup$ Just a usage point: the object you're trying to minimize is a function, not an equation. $\endgroup$ – Matthew Leingang Aug 23 '18 at 19:04
  • $\begingroup$ Edited... Thanks $\endgroup$ – user532063 Aug 23 '18 at 19:06
  • $\begingroup$ I missed the normalized vector constraint. I’ll delete my comment. $\endgroup$ – mathcounterexamples.net Aug 23 '18 at 19:53

Yes, there is a calculus for functions defined in terms of vector and matrix operations. It takes advantage of the algebraic rules for those operations (which is good), but care must be taken because of the noncommutativity of matrix multiplication.

Here's an example. If $x$ is a real variable, $a$ is a constant, and $f(x) = ax^2$, then $f'(x) = 2ax$. if $\mathbf{x}$ is a vector variable, and $A$ is a constant square matrix, then the analogous quadratic function is $f(\mathbf{x}) = \mathbf{x}\cdot A \mathbf{x} = \mathbf{x}^T A \mathbf{x}$. We can express the derivative of $f$ in terms of $\mathbf{x}$.

If the $k$-th component of $\mathbf{x}$ is $x_k$, and the $(i,j)$-entry of $A$ is $a_{ij}$, then $$ f(\mathbf{x}) = \sum_{i,j=1}^n x_i a_{ij} x_j $$ So for each $k$, \begin{align*} \frac{\partial f}{\partial x_k} &= \sum_{i,j=1}^n \frac{\partial}{\partial x_k}(x_i a_{ij} x_j) = \sum_{i,j=1}^n \left(\delta_{ik} a_{ij} x_j + x_i a_{ij} \delta_{jk}\right) \\ &= \sum_{j=1}^n a_{kj}x_j + \sum_{i=1}^n a_{ik}x_i \end{align*} The first summation is equal to the $k$-th component of $A\mathbf{x}$. The second summation is equal to the $k$-th component of $A^T\mathbf{x}$. So we can say: $$ \frac{\partial}{\partial \mathbf{x}}(\mathbf{x}^T A \mathbf{x}) = A \mathbf{x} + A^T \mathbf{x} $$ If $A$ is symmetric, then $\frac{\partial}{\partial \mathbf{x}}(\mathbf{x}^T A \mathbf{x}) = 2A\mathbf{x}$, which is pretty close to $\frac{d}{dx}(ax^2) = 2ax$.

You can read more about this kind of calculus on Wikipedia. The thing is, you almost need to work it out in coordinates first, in order to convince yourself that the matrix algebra shortcuts are valid.

  • 2
    $\begingroup$ So, after inroducing a Lagrange multiplier, you get the equation $2 M_1 v - 4 (\bar{v}^T M_2 v) M_2 v - 2 \lambda v = 0$ as well as $v^T v = 1$. Not so easy to solve... $\endgroup$ – Robert Israel Aug 23 '18 at 19:33
  • $\begingroup$ @RobertIsrael: True $\endgroup$ – Matthew Leingang Aug 23 '18 at 19:35

If you didn't have normalization, this would be exactly solvable in the n=3 case by Sums of Squares methods.

In general this is hard, but there exist approximation algorithms with $1 - Cn^{-2}$ approximation ratios for quartic optimization with quadratic constraints (exactly what you have here).

As far as an explicit form for this problem, I'll write $A=M_1$, $B=M_2$, and $x=v$. The Lagrangian becomes $$L(x,\lambda)=x^* A x - (x^* B x)^2 + \lambda(1 - x^* x)$$ Taking the derivative says that at optimality $$ Ax - 2B x x^* B x - \lambda x=0\text{, or } (A - 2Bxx^*B - \lambda I)x = 0 $$ and, necessarily, $x^*x=1$. As far as I can tell this isn't easily solved.

  • $\begingroup$ Many many thanks, I hope this article really helps... What about matrix differentiation techniques or similar approaches? $\endgroup$ – user532063 Aug 23 '18 at 19:13
  • 1
    $\begingroup$ If you didn't have normalization, the minimimum wouldn't exist: as long as $\bar{v}.M_2.v \ne 0$, you can make the objective an arbitrarily large negative number by appropriate scaling. $\endgroup$ – Robert Israel Aug 23 '18 at 19:20


Note $f(v)=\bar{v}.M_1.v-(\bar{v}.M_2.v)^2$ the function. Your target is to find $v$ such that $f^\prime(v)=0$, where $f^\prime$ is the Fréchet derivative.

To do that you need to use heavily the chain rule writing $f$ as a composition of functions.

For example $\varphi(v)= \bar{v}.M_1v=\varphi_1(\varphi_2(v),v)$ where $\varphi_2(v)=\bar{v}$ and $\varphi_1(u,v)=u.M_1.v$.

Now using Fréchet derivatives,

$$\begin{cases}\varphi_2^\prime(v).h= \bar{h}\\ \varphi_1^\prime(u,v).(h,k)= h.M_1.v+ \bar{u}.M_1k \end{cases}$$ and chain rule

$$\varphi^\prime(v).h=\bar{h}.M_1v +\bar{v}.M_1.h=2\bar{h}.M_1.v$$ using the fact that $M_1$ is hermitian.

You can then follow a similar path for the function $\psi(v)=(\bar{v}.M_2.v)^2$.

And finally use Lagrange multipliers to take into consideration the normalization.

  • $\begingroup$ For the normalization I should add Lagrange multiplier to the equation right? $\endgroup$ – user532063 Aug 23 '18 at 19:43
  • $\begingroup$ Yes indeed. In that case you’ll have to use Lagrange multiplier. $\endgroup$ – mathcounterexamples.net Aug 23 '18 at 19:51
  • $\begingroup$ It is not sufficient that the derivative is zero but also the second derivative must be negative. $\endgroup$ – user532063 Aug 25 '18 at 18:52

Given the matrix $X\in{\mathbb C}^{n\times n},\,\,$ define a function $$F(X) = \begin{bmatrix}{\rm Re}(X) & -{\rm Im}(X)\\{\rm Im}(X) & {\rm Re}(X) \end{bmatrix}\in{\mathbb R}^{2n\times 2n}$$

Similarly, given the vector $x\in{\mathbb C}^{n},\,\,$ define a function $$f(x) = \begin{bmatrix}{\rm Re}(x)\\{\rm Im}(x)\end{bmatrix}\in{\mathbb R}^{2n}$$

Using these functions, your complex problem can be transformed into a real problem, which makes it easier to calculate gradients.

Define some new, real variables $$\eqalign{ A &= F(M_1),\,\,\,B = F(M_2),\,\,\,x = f(v) \cr \alpha &= v^*M_1v = \tfrac{x^TAx}{x^Tx} \cr \beta &= v^*M_2v = \tfrac{x^TBx}{x^Tx} \cr }$$ Since the matrices $(M_1,M_2)$ are hermitian, the matrices $(A,B)$ are symmetric. Also note that $v$ was constrained to be a unit vector, whereas $x$ is unconstrained.

Start by finding the differential and gradient of $\alpha$ $$\eqalign{ d\alpha &= \frac{(2Ax)^Tdx}{x^Tx} - \frac{(x^TAx)(2x^T)dx}{(x^Tx)^2} \cr &= \frac{2(Ax-\alpha x)^Tdx}{x^Tx} \cr \frac{\partial\alpha}{\partial x} &= \tfrac{2(A-\alpha I)x}{x^Tx} \cr }$$ Similarly, the gradient of $\beta$ is $\frac{\partial\beta}{\partial x} = \frac{2(B-\beta I)x}{x^Tx}$

The gradient of your target function can be calculated as $$\eqalign{ \phi &= \alpha-\beta^2 \cr \frac{\partial\phi}{\partial x} &= \frac{\partial\alpha}{\partial x} - 2\beta\frac{\partial\beta}{\partial x} = \frac{2\,\big(Ax-2\beta Bx + 2\beta^2x-\alpha x\big)}{x^Tx} \cr }$$ Now use your favorite gradient-based unconstrained minimization method (e.g. Barzilai-Borwein) to solve for $x$ -- and remember that $(\alpha,\beta)$ are not constants, but are functions of $x$.

Once you have $x$, you can recover the $v$ vector which solves the original complex problem.

Here's some Julia code

function q(A,x); return sum(x'*A*x); end        # quadratic form
function normF(x); return sqrt(sum(x'*x)); end  # frobenius norm

# objective function and gradient
function     F(A,B,x); return q(A,x) - q(B,x)^2; end
function gradF(A,B,x)
  a,b,c = q(A,x), q(B,x), q(one(A),x)
  return (2/c) * ((A*x-a*x) - 2*b*(B*x-b*x))

# Barzilai-Borwein
function bbSolve(tol,nMax,A,B,x)
   # initialize
   g=x.+0; dx=x*0;  dg=x*0; xBest=x.+0; g=gradF(A,B,x); vBest=v=normF(g)
   n=0; b=1e-5; dx=b*(x.+b); dg=g.+0; g=gradF(A,B,x-dx); dg -= g

   # iterate
   while true
       n += 1
       b = sum(dg'*dx) / sum(dg'*dg)  # barzilai steplength
       dx = -b*g; x += dx
       dg = -g; g=gradF(A,B,x); dg += g
       v = normF(g)

       if v < vBest
          vBest = v; xBest = x

       if isnan(v) || isinf(v)
       elseif v < tol || n > nMax
   return n, xBest

# test with randomly generated matrices
n=5;    x=ones(2*n,1);     # initial guess x = 1
M1=rand(n,n)+im*rand(n,n); M1+=M1';  # hermitian
M2=rand(n,n)+im*rand(n,n); M2+=M2';
A=[real(M1) -imag(M1); imag(M1) real(M1)];
B=[real(M2) -imag(M2); imag(M2) real(M2)];
@time k,x = bbSolve(1e-14, 200, A,B,x);
@printf( "k, F(x), |x|, x  =  %s\n",  (k, F(A,B,x), normF(x), x) );

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy