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](https://mathoverflow.net/questions/29958/eigenvalue-problems-with-linear-constraints) is the same, but I couldn't parse the actual question to verify duplicity).



Linearly constrained eigenvalue problem
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.


