Analytical solution for equality-constrained QCLP

Question

I'm looking for a way to derive a closed form expression for $\omega$, which is the solution of the following maximization problem: $$\max_{w} \omega^{T}\mu$$ $$\textrm{st.}\hspace{0.5cm} \omega^{T}\Sigma\omega=\sigma^2$$ $$ \omega^{T}\textbf{1}=1 $$ In particular I'm struggeling at the point where i try to solve for the lagrange multipliers given the second FOC: $$\frac{d{L(\dots)}}{d{\omega}}=\mu-2\lambda_{1}\Sigma\omega-\lambda_{2}\textbf{1}=0$$ $$\frac{d{L(\dots)}}{d{\lambda_{1}}}=\sigma^{2}-\omega^{T}\Sigma\omega=0$$ $$\frac{d{L(\dots)}}{d{\lambda_{2}}}=1-\omega^{T}\textbf{1}=0$$

Here, $\mu$ is a (1xn)-vector nx1 of some constants, $\omega$ the (1xn)-vector of interest, $\Sigma$ some (nxn)-matrix, $\lambda_{1}$ and $\lambda_{2}$ the lagrange multipliers for constraints 1 and 2 and $\sigma^2$ is some (1x1)-variable.

The problem here (at least for me) is that the first derivative for the lagrange multipliers of the constraints includes two $\omega$-vectors and thus yields a quadratic expression, for which i can't easily solve the second FOC for the needed lagrange multipliers.. I'm not sure whether i'm on the right track and if yes, how to solve for the analytical expression of $\omega$ from here.

If an analytical solution / expression for $\omega$ does exist, how can I explicitly derive this solution step-by-step?

Answer 1

Here's another solution that doesn't reduce the order of the problem.

We have the following three equations (as described in the question)

$ \mu - 2 \lambda_1 \Sigma \omega - \lambda_2 \mathbf{1} = \mathbf{0} $

$ \omega^T \Sigma \omega = \sigma^2 $

$ \omega^T \mathbf{1} = 1 $

From the first equation,

$ \omega = \dfrac{1}{2 \lambda_1} \Sigma^{-1} ( \mu - \lambda_2 \mathbf{1} ) $

substituting this into the next two equations, gives us,

$\dfrac{1}{(2 \lambda_1)^2} ( \mu - \lambda_2 \mathbf{1} )^T \Sigma^{-1} ( \mu - \lambda_2 \mathbf{1} ) = \sigma^2 $

and

$\dfrac{1}{2 \lambda_1} \mathbf{1}^T \Sigma^{-1} ( \mu - \lambda_2 \mathbf{1} ) = 1 $

Combining the two equations by eliminating $(2 \lambda_1)$,

$ ( \mu - \lambda_2 \mathbf{1} )^T \Sigma^{-1} ( \mu - \lambda_2 \mathbf{1} ) = \sigma^2 ( \mathbf{1}^T \Sigma^{-1} ( \mu - \lambda_2 \mathbf{1} ) )^2 $

This is a quadratic equation in $\lambda_2$ that can be solved for $2$ values of $\lambda_2$. Corresponding to each value of $\lambda_2$ we then solve for $\lambda_1$, and $\omega$ follows.

I will now expand the above quadratic equation and put it in the standard quadratic equation format.

$ \lambda_2^2 (\mathbf{1}^T \Sigma^{-1} \mathbf{1} ) - 2 \lambda_2 (\mu^T \Sigma^{-1} \mathbf{1} ) + \mu^T \Sigma^{-1} \mu = \sigma^2 ( \mu^T \Sigma^{-1} \mathbf{1} - \lambda_2 (\mathbf{1}^T \Sigma^{-1} \mathbf{1} ) )^2 $

Define,

$ a = \mathbf{1}^T \Sigma^{-1} \mathbf{1} $

$ b = \mu^T \Sigma^{-1} \mathbf{1} $

$ c = \mu^T \Sigma^{-1} \mu $

Then our equation in $\lambda_2$ becomes

$ a \lambda_2^2 - 2 b \lambda_2 + c = \sigma^2 ( b - a \lambda_2 )^2 $

Expanding the squared term on the right hand side and re-grouping,

$ (a - a^2 \sigma^2 ) \lambda_2^2 + ( 2 a b \sigma^2 - 2 b ) \lambda_2 + (c - b^2 \sigma^2) = 0 $

Divide through by $ (1 - a \sigma^2 ) $,

$ a \lambda_2^2 - 2 b \lambda_2 + \bigg(\dfrac{ c - b^2 \sigma^2}{1- a \sigma^2 } \bigg) = 0 $

Complete the square,

$ a ( \lambda - \dfrac{b}{a} )^2 + \bigg(\dfrac{ c - b^2 \sigma^2}{1- a \sigma^2 } \bigg) - \dfrac{b^2}{a} = 0 $

The constant term simplifies to

$ \dfrac{ a (c - b^2 \sigma^2 ) - b^2 (1 - a \sigma^2 ) }{a (1 - a \sigma^2)} = \dfrac{ ac - b^2 }{ a (1 - a \sigma^2) } $

Divide through by $a$, then

$ ( \lambda - \dfrac{b}{a} )^2 = \dfrac{ b^2 - a c }{ a^2 (1 - a \sigma^2)} $

Hence, the solutions are

$ \lambda_2 = \dfrac{b}{a} \pm \sqrt{ \dfrac{ b^2 - a c }{a^2 (1 - a \sigma^2) } } $

Note that $a$ and $c$ are positive. To have a real solution, we require that

$ \dfrac{ b^2 - a c }{a^2 (1 - a \sigma^2) } \ge 0 $

Now, for the tangent plane to the ellipsoid $ \omega^T \Sigma \omega = \sigma^2 $ to have a normal vector along $\mathbf{1}$, we have

$ \Sigma \omega = \alpha \mathbf{1} $

Substituting in the ellipsoid equation, we get

$ \alpha = \sqrt{\dfrac{ \sigma^2 }{ \mathbf{1}^T \Sigma^{-1} \mathbf{1} } } = \sqrt{\dfrac{ \sigma^2 }{a} } $

So that tangency point is

$ \omega_1 = \dfrac{ \Sigma^{-1} \mathbf{1} \sqrt{ \sigma^2} }{ \sqrt{a} } $

The distance between the origin and the tangent plane is

$ d_1 = \sqrt{\sigma^2} \dfrac{ \mathbf{1}^T \Sigma^{-1} \mathbf{1} } { \sqrt{a} \sqrt{\mathbf{1}^T \mathbf{1} } } = \dfrac{ \sqrt{a \sigma^2} }{\sqrt{n}} $

On the other hand, the distance between the origin and the hyperplane $ \omega^T \mathbf{1} = 1 $ is given by

$ d_2 = \dfrac{1}{\sqrt{ n } } $

Hence, there will be an intersection between the ellipsoid and the hyperplane if and only if $d_1 \gt d_2 $, which is equivalent to

$ a \sigma^2 \gt 1 $

Using this in the inequality for real solutions of the quadratic in $\lambda_2$ above, we deduce that the condition is now

$ b^2 - a c \le 0 $

i.e

$ a c \ge b^2 $

Careful examination of the definitions of $a,b, c$, and the use of the Cauchy-Schwartz inequality, leads to the conclusion that we always have,

$ a c \ge b^2 $

Therefore, in conclusion, the necessary and sufficient condition that the quadratic equation has real solutions is that

$ a \sigma^2 \ge 1 $

And this inequality is satisfied if and only if the hyperplane $\omega^T \mathbf{1} = 1$ intersects the ellipsoid $\omega^T \Sigma \omega = \sigma^2$.

The two values of $\lambda_2 $ are given by

$ \lambda_2 = \dfrac{1}{a} \bigg( b \pm \sqrt{ \dfrac{ ac - b^2 }{ a \sigma^2 - 1 } } \bigg) $

One solution is for the minimum and the other the maximum of the given objective function.

From the value of $\lambda_2$ (there is obviously two of them), we can calculate the corresponding $\lambda_1$ from

$\dfrac{1}{2 \lambda_1} \mathbf{1}^T \Sigma^{-1} ( \mu - \lambda_2 \mathbf{1} ) = 1 $

which is stated above in the derivation, and from which, we have

$ \lambda_1 = \dfrac{1}{2} \mathbf{1}^T \Sigma^{-1} (\mu - \lambda_2 \mathbf{1} ) = \dfrac{1}{2} ( b - a \lambda_2 ) $

And finally, from $\lambda_1, \lambda_2 $ we can compute $\omega$.

This is given by,

$ \omega = \dfrac{1}{2 \lambda_1} \Sigma^{-1} ( \mu - \lambda_2 \mathbf{1} ) $

Now

$ f(\omega) = \mu^T \omega = \dfrac{1}{2 \lambda_1} ( c - \lambda_2 b ) $

Substituting the above expression for $\lambda_1$, we get

$ f = \dfrac{ c - b \lambda_2 }{ b - a \lambda_2 } $

Answer 2

First reduce the order of the system by 1, by solving the linear system

$ \omega^T \mathbf{1} = 1 $

This will give

$ \omega = \omega_0 + M u $

where $u \in \mathbb{R}^{n-1} $ and $ M \in \mathbb{R}^{n \times (n-1) } $

Now the problem is re-written as follows

$$ \displaystyle \max\limits_{u} \mu^T M u $$

Subject to

$$ (\omega_0 + M u)^T \Sigma (\omega_0 + M u) = \sigma^2 $$

by defining $ \mu_1 = M^T \mu $ and $ A = M^T \Sigma M, \\B = 2 M^T \Sigma \omega_0,\\ C = \omega_0^T \Sigma \omega_0 - \sigma^2 \\ \\$

the problem becomes

$$ \displaystyle \max\limits_{u} \mu_1^T u $$

Subject to

$$ u^T A u + u^T B + C = 0 $$

The gradient of the Lagrangian is

$ \nabla_u = \mathbb{0} = \mu_1 + \lambda (2 A u + B ) $

From here, we get,

$ u = \dfrac{1}{2} A^{-1} ( - B - \dfrac{1}{\lambda} \mu_1 ) $

Plug this into the constraint equation,

$ u^T A u + u^T B + C = 0 $

Take $t = \dfrac{1}{\lambda} $, then

$ \dfrac{1}{4} ( - B - t \mu_1 )^T A^{-1} (- B - t \mu_1 ) + \dfrac{1}{2} B^T A^{-1} (- B - t \mu_1 ) + C = 0 $

This is a quadratic equation in $t$, which can be solved readily by the quadratic formula.

From the two roots $t$, we have two values for $\lambda = \dfrac{1}{t} $ , and $u$ follows, and so does $\omega$. The two values of $t$ correspond to the minimum and maximum values of the function.

Here's a subroutine written in Microsoft VB for Applications, to solve $a^T x = b$ given vector $a$ and scalar $b $.

  Public Sub solve_hyperplane_equation()
  Dim a(10), b As Double
  Dim x0(10), M(10, 10) As Double
  Dim i, j, k As Integer

  n = 5

  a(1) = 1
  a(2) = -2
  a(3) = 0.5
  a(4) = 0
  a(5) = 4

  b = 10

  ' we want to solve a^T x = b

  ' The solution is of the form x = x0 + M u

  ' where the columns of V are orthogonal to a

  ' first find x0

  ' augment vector a with the scalar b

  a(n + 1) = b

  pivotfound = 0

  For i = 1 To n
      If Abs(a(i)) > 0.00000001 Then
         pivotfound = 1
         t = a(i)
         i0 = i
         Exit For
      End If
  Next i

  If pivotfound = 0 Then
  MsgBox ("Invalid input")
  End If

  t = a(i0)

  For j = 1 To n + 1
  a(j) = a(j) / t
  Next j

  ' the pivot is i0, all the other variables are arbitrary scalars

  ' find x0

  ' set all varibles to zero , except x_i0, solve for it

  For k = 1 To n
  x0(k) = 0
  Next k

  x0(i0) = a(n + 1)

  icount = 0

  For k = 1 To n

  If k = i0 Then GoTo lnext_k

  icount = icount + 1

  ' a(1) x1 + a(2) x2 + .... + a(i0) x_i0 + .... + a(n) x_n = a(n+1)

  ' set all the arbitrary variables to zero, except variable k, set to 1, solve for x_i0

  M(k, icount) = 1

  M(i0, icount) = (-a(k)) / a(i0)

  lnext_k:
  Next k

  ' display results

  For k = 0 To n - 1

  For i = 1 To n

  If k = 0 Then ActiveSheet.Cells(i, k + 1) = x0(i) Else ActiveSheet.Cells(i, k + 2) = xx(i, k)

  Next i
  Next k

  For i = 1 To n + 1
  ActiveSheet.Cells(n + 2, i) = a(i)
  Next i

  End Sub

The output of the above code follows