# How can I numerically solve a saddle point problem with repeated constraints?

I am interested in numerically solving the following constrained minimization problem; Find the value of $x\in \mathbb{R}^n$ that minimizes $f$ where $f\colon \mathbb{R}^n\to \mathbb{R}$ is defined as $$f(x) = \frac{1}{2}(x,A x)-(x,b),$$ for the matrix $A$ being large sparse, symmetric and positive definite, subject to $$B^\mathrm{T} x = g.$$

Here the matrix $B\in \mathrm{R}^{n\times m}$ is large sparse, and $m$ is large but small compared to $n$. So far so good, this is just the normal saddle point problem. The wrinkle here is that $B$ may have linearly dependent columns, and that I know $g\in \mathrm{Range}(B^\mathrm{T})$, at least up to machine precision.

The problem comes from a finite element code where I generate the constraint matrices by against too many test functions. I can explain more if anyone is interested. (I Should be testing against all elements in trace of some space, S, but I don't have that space available to me easily so I test against all elements in S.

How can I efficiently solve this problem? I am mainly interested in some sort of Krylov / Relaxation method. Below I outline a method using the SVD of $B$, and a $QR$ method would work too.

This is what I would do for small matrices: (Spoiler Alert: this is very tedious and somewhat boring, I wrote it out in this much detail for my own reference, and to show that I did some work before asking here.)

Let $(U\in \mathbb{R}^{n\times n},\Sigma\in \mathbb{R}^{n\times m}, V\in \mathbb{R}^{m\times m})$ constitute the singular value decomposition of $B$, so that $$U\Sigma V^\mathrm{T} = B.$$

If $r$ is the rank of $B$, then the first $r$ columns of $U$ (call then $U_1$) span the range of $B$ and the last $n-r$ columns of $U$ (call them $U_2$) span the nullspace of $B^\mathrm{T}$. Similarly, the first $r$ columns of $V$ (call then $V_1$) span the range of $B^\mathrm{T}$ and the last $m-r$ columns (call them $V_2$) span the null space of $B$. We also define $\Sigma_1$ as $r\times r$ matrix that sits at tthe top-left hand corner of $\Sigma$.

We now decompose $x$ in $x = \hat{x} + x_\perp$ where $\hat{x} \in \mathrm{Null}(B^\mathrm{T})$ and $x_\perp\in (\mathrm{Null}(B^\mathrm{T}))^{\perp}= \mathrm{Range}(B)$. This means that $x_\perp = U_1\alpha$ for some $\alpha\in \mathrm{R}^{n-r}$ and that $\hat{x} = U_2\beta$ for some $\beta\in \mathbb{R}^{n-t}$.

Subbing this second identity into the constraint we get $$V\Sigma^\mathrm{T} U^\mathrm{T} U_1\alpha = g,$$

We can now multiply on the left by $V_1^\mathrm{T}$ which gives $V_1^\mathrm{T}V\Sigma^\mathrm{T} U^\mathrm{T} U_1\alpha = V_1^\mathrm{T}g$, and noting that $V_1^\mathrm{T}V\Sigma^\mathrm{T} U^\mathrm{T} U_1 = \Sigma_1^\mathrm{T}=\Sigma_1$ we get $\alpha = \Sigma_1^{-1} V_1^{\mathrm{T}} g$ and $x_\perp$ is known.

Note that if we had multiplied on the left by $V_2^\mathrm{T}$ we have obtained $0 = 0$ since $g\in \mathrm{Range}(B^\mathrm{T})$ so all the "information" in the constraints was retained when we multiplied on the left by $V_1^\mathrm{T}$.

Then we define $F:\mathrm{R}^{n-r}\to \mathrm{R}$ as $$F(\beta) = \frac{1}{2}(x_\perp + U_2\beta,A(x_\perp+U_2\beta) ) - (b, x_\perp + U_2\beta)$$

Now we find the $\beta$ that minimizes this functional (symmetric positive definite) and we can find $\hat{x}$.

• Which method would you use for the saddle point matrix if matrix $B$ is full rank? – VorKir Jul 3 '17 at 19:49