If I have the problem $$\min_{x,y} \|Ax - b\|_2^2 \quad s.t. \quad Cx = D^T y$$ (for $A\in \mathbb{R}^{m\times n}$ with $m < n$, $C \in \mathbb{S}^{n}$, $D \in \mathbb{R}^{k \times n}$ full row rank $k < n$ orthogonal projection, $x \in \mathbb{R}^n$, $y \in \mathbb{R}^k$, and $b \in \mathbb{R}^m$), how can I formulate the problem as a system of linear equations to solve? I want to enforce a the relationship between $x$ and $y$ as a projection constraint but retain both $x$ and $y$ in the formulation (i.e., not eliminate $y$ by requiring the orthogonal projection $(I- D^T D) Cx = 0$).
My Lagrangian $\mathcal{L}(x,y,\lambda) = x^T A^T Ax - 2x^T A^Tb + \lambda^T(Cx - D^T y)$ gives:
$$ \begin{align} \nabla_x \mathcal{L} = 0 &\implies A^T Ax + C^T \lambda = A^T b \\ \nabla_y \mathcal{L} = 0 &\implies D\lambda = 0 \\ \nabla_x \mathcal{L} = 0 &\implies Cx = D^Ty \\ \end{align} $$
I tried to put this into a single, symmetric system so that I could use MINRES:
$$ \begin{bmatrix} A^T A & C^T & 0 \\ C & -D^T & 0 \\ 0 & 0 & D \end{bmatrix} \begin{bmatrix} x \\ y \\ \lambda \end{bmatrix} = \begin{bmatrix} A^T b \\ 0 \\ 0 \\ \end{bmatrix} $$ but my $D^T$ dimension doesn't fit within the dimensions of $C$ and $C^T$. What am I missing / doing wrong?
We have the following constrained least-squares problem in $\mathrm x_1 \in \mathbb R^{n_1}$ and $\mathrm x_2 \in \mathbb R^{n_2}$
$$\begin{array}{ll} \text{minimize} & \| \mathrm A_1 \mathrm x_1 - \mathrm b \|_2^2\\ \text{subject to} & \mathrm C_1 \mathrm x_1 + \mathrm C_2 \mathrm x_2 = \mathrm 0\end{array}$$
We rewrite the original problem as follows
$$\begin{array}{ll} \text{minimize} & \| \mathrm A \mathrm x - \mathrm b \|_2^2\\ \text{subject to} & \mathrm C \mathrm x = \mathrm 0\end{array}$$
where $\mathrm x := \begin{bmatrix} \mathrm x_1\\ \mathrm x_2\end{bmatrix}$ and $\mathrm A := \begin{bmatrix} \mathrm A_1 & \mathrm O\end{bmatrix}$ and $\mathrm C := \begin{bmatrix} \mathrm C_1 & \mathrm C_2\end{bmatrix}$. Let the Lagrangian be
$$\mathcal L (\mathrm x, \mu) := \frac 12 \| \mathrm A \mathrm x - \mathrm b \|_2^2 + \mu^{\top} \mathrm C \mathrm x$$
Taking the partial derivatives and finding where they vanish, we obtain the following linear system
$$\begin{bmatrix} \mathrm A^\top \mathrm A & \mathrm C^\top\\ \mathrm C & \mathrm O\end{bmatrix} \begin{bmatrix} \mathrm x\\ \mu\end{bmatrix} = \begin{bmatrix} \mathrm A^\top \mathrm b\\ 0\end{bmatrix}$$
which can be rewritten as follows
$$\begin{bmatrix} \mathrm A_1^\top \mathrm A_1 & \mathrm O & \mathrm C_1^\top\\ \mathrm O & \mathrm O & \mathrm C_2^\top\\ \mathrm C_1 & \mathrm C_2 & \mathrm O\end{bmatrix} \begin{bmatrix} \mathrm x_1\\ \mathrm x_2\\ \mu\end{bmatrix} = \begin{bmatrix} \mathrm A_1^\top \mathrm b\\ 0\\ 0\end{bmatrix}$$