Linearly constrained eigenvalue problem

Question

Suppose I'd like to:

\begin{align} \mathop{\text{min}}_\mathbf{x} && \mathbf{x}^T\mathbf{A}\mathbf{x} \\ \text{subject to:} && \mathbf{x}^T \mathbf{M} \mathbf{x} = 1\\ && \mathbf{C}\mathbf{x} = \mathbf{b} \end{align} where all vector variables are known except $\mathbf{x}$, and $\mathbf{C}$ is full row rank.

If I didn't have the $\mathbf{C}\mathbf{x} = \mathbf{b}$ constraint then after applying the Lagrange multiplier method I could solve this as a generalized eigenvalue problem:

\begin{align} \text{solve} && \mathbf{A}\mathbf{x} = \lambda \mathbf{M} \mathbf{x}\\ \text{subject to:} && \mathbf{x}^T \mathbf{M} \mathbf{x} = 1\\ \end{align}

If I try to apply the same approach to my original problem I get something that doesn't look quite like a generalized Eigenvalue problem:

\begin{align} \text{solve} && \left(\begin{matrix} \mathbf{A} & \mathbf{C}^T\\ \mathbf{C} & \mathbf{0} \end{matrix}\right) \left(\begin{matrix} \mathbf{x} &\\ \mu \end{matrix}\right) = \left(\begin{matrix} \lambda \mathbf{M} \mathbf{x}\\ \mathbf{b} \end{matrix}\right)\\ \text{subject to:} && \mathbf{x}^T \mathbf{M} \mathbf{x} = 1 \\ \end{align}

Ideally, I'd like to reduce my problem to an instance of the generalized eigenvalue problem so I can use an off-the-shelf numerical solver.

What's the best way to solve this problem?

(The title of this question is the same, but I couldn't parse the actual question to verify duplicity).

Answer 1

I think Pushpendre's answer isn't quite right, but it gets you most of the way there. Getting rid of that pesky constant term is a bit tricky relative to the homogeneous case.

Let's take his suggested substitutions: $$ \begin{array}{rl} M&:=N^\top N\\ y&:=Nx\\ D&:=CN^{-1}\\ B&:=N^{-\top}AN^{-1} \end{array} $$ I don't think these are strictly necessary (e.g. this encodes an assumption that $M$ is semidefinite), but it simplifies notation quite a bit for the most realistic/common case. I'm going to be cavalier about assuming certain matrices are symmetric/invertible as convenient.

Then, your problem becomes $$ \begin{array}{rl} \min_y &y^\top By\\ \mathrm{s.t.} & \|y\|^2=1\\ & Dy=b \end{array} $$ This problem has Lagrange multiplier expression (I am sprinkling in constant factors so that they simplify my arithmetic later) $$ \Lambda(y;\lambda,\mu):=\frac{1}{2}y^\top By+\frac{1}{2}\lambda(1-\|y\|^2)+\mu^\top (b-Dy) $$ Differentiating with respect to $y$ shows $$ 0=\nabla_y\Lambda(y;\lambda,\mu)=By-\lambda y-D^\top\mu. $$ Here I am assuming w.l.o.g. that $B$ is symmetric.

Pre-multiplying the critical point condition by $D$ shows $$ DBy=\lambda Dy + DD^\top\mu=\lambda b + DD^\top \mu. $$ Let's further assume that $D$ has full rank. The most common case is that $D\in\mathbb{R}^{m\times n},$ where $m<n$; I'll assume $D$ has rank $m$. Then, $DD^\top$ is an invertible matrix, and $D$ admits a pseudoinverse $D^+:=D^\top (DD^\top)^{-1}$ satisfying $DD^+=I$.

Then, we can isolate $\mu$ as $$ \mu=(DD^\top)^{-1}(DBy-\lambda b). $$ Plugging this back into the expression for $\nabla_y\Lambda$ shows $$ \begin{array}{rl} 0&=By-\lambda y-D^\top (DD^\top)^{-1}(DBy-\lambda b)\\ &=By-\lambda y-D^+(DBy-\lambda b)\\ \implies [(I-D^+D)B-\lambda I]y&=-\lambda D^+b. \end{array} $$ For each $\lambda$, this expression gives a system of equation solvable for $y$. In other words, we can think of $y$ as a function $y(\lambda)$ of $\lambda$.

Define a function $f(\lambda):\mathbb{R}\rightarrow\mathbb{R}\cup\{-\infty,\infty\}$ as $$f(\lambda):=\|y(\lambda)\|^2-1.$$ Use any standard 1d method to find a root $\lambda$ such that $f(\lambda)=0.$ Then, $y(\lambda)$ is a critical point of your optimization problem. To be totally rigorous, you should check that $y(\lambda)$ satisfies the constraint $Dy=b$, but this is easy to check for any $\lambda$ from the system for $y(\lambda)$ (multiply both sides by $D$).

Note that numerically the problem of finding roots of $f(\cdot)$ isn't great. $f$ likely has multiple roots and asymptotes where the system of equations for $y$ is not solvable. It is known as a secular equation, for which some specialized solution algorithms exist.

---

NOTE: This trick is similar to the one suggested for "LSQI" in Golub/Van Loan's book.

Answer 2

Your problem has been answered at https://scicomp.stackexchange.com/questions/14096/sparse-smallest-eigenvalue-problem-on-a-linear-subspace :) Or you can read Golub's original paper Some modified matrix eigenvalue problems

The basic intuition is that basically you want to find the eigenvalues over a subspace that is defined $Cx = b$ so just find its basis and convert the rayleigh quotient to that basis.

EDIT: I just realized that you need the space which satisfies $Cx=b$ and zero need not be in this space so its not a subspace. But I think the techniques in those answers could still be adapted.

WRONG One way I can think of handling the nonzero right hand size is to make an augmented matrix $C'=[C\; b]$ and augmented $x' = [x; 1]$ then $C'x' = 0$ and then also create $M' = [M\; 0; 0\; 0]$ and same for $A$. After this augmentation and zero padding everything works out.

2nd try: we can re write the objective as $$\arg \min \frac{X^TAX}{X^TMx} \textrm{subject to } C(x-x')=0$$ where $x'$ is any point satisfying $Cx=b$ and then assuming that $M$ is positve definite which means that it can be broken into $M = N^TN$ (assuming real.)

let $y = Nx$.

let $B = (N^{-1})^TAN^{-1}$.

let $D = CN^{-1}$.

then the objective becomes $$\arg\min \frac{y^TBy}{y^Ty} \textrm{ subject to } D(y-y') = 0$$. Now we are pretty close to Golub's setup in 1.1, 1.2 and 1.3 but not quite because of the non-zero RHS. But we can still use the lagrangian (assuming everywhere things were positive definite and inverses were well conditioned (practically))

So the lagrangian is $$y^TBy - \lambda (y^Ty -1) + 2 \mu^T(D(y-y'))$$

EDIT: See Justin's answer above for a correct explanation

The first derivative is $$By - \lambda y + 2 \mu^TD = 0$$ multiply with $D^T$ on the right and use the constraint that $D^Ty = D^Ty'$ to get $$By'B^T - \lambda y'B^T + 2 \mu^T DD^T = 0$$ which gives us a value of $\mu$ (because we know $B$ and $y'$ and I assume that since $C$ was full row rank therefore $DD^T$ would be invertible, if not invertible then we have a range of solutions for $\mu$ and would have to pick the best) then substitute to get a generalized eigen value problem in terms of lambda which would be the solution.