I am solving a system of the form $$ \begin{pmatrix} A & b^T \\ b & 0 \end{pmatrix} \begin{pmatrix} x \\ \ell \end{pmatrix} = \begin{pmatrix} c\\ 0 \end{pmatrix} $$ Where $A$ is a symmetric matrix with a one dimensional kernel. $A$ is positive definite on the orthogonal of the kernel. $b$ is a vector and $c$ is such that $Ax=c$ has solutions.

The matrix $A$ comes from the finite element discretization of $$ -\Delta U - \lambda_1 U = f $$ with $U=0$ on the boundary of the domain and $\lambda_1$ is the first eigenvalue of the Dirichlet Laplacian. Of course, the previous equation has infinitely many solutions: if $U$ is a solution, so is $U+\alpha u_1$ ($u_1$ is the eigenfunction associated to $\lambda_1$). Therefore adding an orthogonality constraint on the first eigenfunction $u_1$ is natural: $\int_\Omega Uu_1 = 0$.

If $K,M$ are the rigidity and mass matrices, I have $A = K-\lambda M$, $B = Mu_1$. The right hand side verifies a compatibility condition $(f,u_1) = 0$ (theoretically and discretely). Therefore, I know that the Lagrange multiplier $\ell$ is zero.

I have an efficient version for solving the above system using Petsc and fieldsplit preconditioner. I built a smaller test case in Matlab and I noticed that using the regular backslash is slow compared with a naive Schur complement I implemented. Therefore my question:

Is there a standard way to reduce the above system to a simpler one for which standard solvers work well?

Are there other efficient ways of solving $Ax=c$ when $\dim \ker A = 1$ and an arbitrary one dimensional orthogonality constraint is considered? (and when we know that solutions exist)

Thank you in advance!


1 Answer 1



You mention this:

Are there other efficient ways of solving $Ax=c$ when $\dim\ker A=1$ and an arbitrary one dimensional orthogonality constraint is considered? (and when we know that solutions exist)

But you also considered a case that results from a minimisation problem, so I'll go through both.

System Constrained Case

Let $u_1$ be the only eigenvector of $A$ corresponding to a zero eigenvalue. Then the strong form of your constraint reads:

$$\begin{bmatrix} A \\ u_1^T \end{bmatrix} x = \begin{bmatrix} c \\ 0 \end{bmatrix}.$$

If $x$ satisfies $u_1^T x = 0$ then it is from the kernel of $u_1^T$. A projection onto the kernel of $u_1^T$ is given as $Pr_{\ker(u_1^T)} = (I-u_1u_1^T)$ (if the vector was not normalized you would have used $I-\frac{u_1u_1^T}{u_1^Tu_1}$). This means that you can just solve $A(I-u_1u_1^T)x = c$. However, because $u_1$ is in the kernel of $A$ then $A(I-u_1u_1^T)x = Ax = c$. So technically you don't need to do anything as long as you have a solver that can handle a singular system. The conjugate gradient can do so if $A$ is symmetric positive semi-definite. In fact if you set the initial guess $x_0$ to be such that $u_1^Tx_0=0$ (for example if you set it to zero) then you will get the desired Moore-Penrose pseudoinverse solution. If you set it to a vector that has a non-zero part in the kernel, then it is sufficient to take $(I-u_1u_1^T)x$ instead.

Note that the above will have no solution if $c$ is not in the range of $A$. This is actually a difference to the minimisation problem in the next section, which always has a solution.

Let's consider the more general problem: $$\begin{bmatrix} A \\ u_1^T \end{bmatrix} x = \begin{bmatrix} c \\ \mu \end{bmatrix}.$$

Note that $u^T_1x = \mu$ is equivalent to $u_1u_1^Tx = u_1 \mu$. At the same time you have $Ax = c$. Since $u_1$ is not in the range of $A$, and $c$ is assumed to not be in the kernel of $A$ you can just add these two equations without modifying the set of solutions $Ax + u_1u_1^Tx = c+u_1\mu$. But $\tilde{A} = A+u_1u_1^T$ is a full-rank matrix in your case. This means you could also interpret this as just solving the problem $\tilde{A}x = c+u_1\mu$. In practice this results in the same with the conjugate gradient solver as solving $Ax = c$ and just setting the initial guess $x_0 = u_1\mu$. So you really don't have to deal with any constraints.

Minimisation of the Constrained Case

Alternatively one could formulate the minimisation problem: $$\min_{u_1^Tx = 0}\frac{1}{2}x^TAx - x^Tc.$$ Rewriting it as an unconstrained problem with Lagrange multipliers yields: $$\max_{\ell} \min_x \frac{1}{2}x^TAx - x^Tc + \ell u_1^T x.$$ This results in the system you wrote in the beginning: $$\begin{bmatrix} A & u_1 \\ u_1^T & 0 \end{bmatrix} \begin{bmatrix} x \\ \ell \end{bmatrix} = \begin{bmatrix} c \\ 0 \end{bmatrix}.$$

You can however impose the constraint earlier: $$\min_{u_1^Tx = 0}\frac{1}{2}x^TAx - x^Tc \implies \min_x \frac{1}{2}x^T(I-u_1u_1^T)A(I-u_1u_1^T)x - x^T(I-u_1u_1^T)c.$$ This results in the system: $$(I-u_1u_1^T)A(I-u_1u_1^T)x = (I-u_1u_1^T)c \implies Ax = (I-u_1u_1^T)c.$$ Unlike the system constrained case this always has a solution since it guarantees that $(I-u_1u_1^T)c$ is in the range of $A$ (at least as long as $u_1$ is the only eigenvector of $A$ that corresponds to a zero eigenvalue, since then $I-u_1u_1^T$ is a projection onto the image of $A$). But if $c$ was in $\img(A)$ to begin with, then of course $c = (I-u_1u_1^T)c$, so the two are equivalent in such a case.

System General Constrained Case

I also derived the general case just in case you need it.

Consider the problem $Ax=a$ under the constraint that $Bx=b$. Suppose $B$ has the singular value decomposition: $$B = \begin{bmatrix} P_{\img} & P_{\ker} \end{bmatrix}\begin{bmatrix} D & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} Q_{\img}^* \\ Q_{\ker}^* \end{bmatrix}.$$ Set $x = Q_{\img}z_{\img} + Q_{\ker} z_{\ker}$ and plug it in $Bx=b$, then: $$Bx = P_{\img}Dz_{\img} = b \implies z_{\img} = D^{-1}P_{\img}^*b.$$ Note that the implication is valid only if $Bx=b$ has a solution. Otherwise the above "solution" solves the projected problem $Bx = b^+ = P_{\img}P^*_{\img}b$, where we have projected $b$ onto the image of $B$. In either case, term $B^+_{\img} = Q_{\img}D^{-1}P_{\img}^*$ then $x_B = B^{+}_{\img}b = Q_{\img} z_{\ker}$. We can now plug $x = x_B + Q_{\ker} z_{\ker}$ into $Ax=a$ and solve for $z_{\ker}$: $$AQ_{\ker}z_{\ker} = a-Ax_B.$$ Note that $Q_{\ker}z_{\ker}$ is a vector in the kernel of $B$. If you want to directly work with $x$, we can instead solve the following system $APr_{\ker(B)}x = a-Ax_B$, where $Pr_S$ is a projection onto $S$. This can be implemented as $Pr_{\ker(B)} = Q_{\ker}Q_{\ker}^* = I-Q_{\img}Q_{\img}^*$ if $Q_{\ker}$ or $Q_{\img}$ are readily available. If they are not, one could use $Pr_{\ker(B)} = I-B^+B$ where the application of the Moore-Penrose pseudoinverse $B^+$ could be computed e.g. using the conjugate solver for the normal equations with an initial guess of zero (since that converges to the Moore-Penrose pseudoinverse solution).

Also in relation to your problem, note that if $Q_{\img}$ is in the kernel of $A$ then of course $Ax_B = 0$ and $APr_{\ker(B)}x = Ax$, which recovers what I wrote before as a special case.

Minimisation of the General Constrained Case

Consider the more general constraint: $$\min_{Bx=b} \frac{1}{2}x^TLx - x^Ta.$$

Suppose you have the singular value decomposition of $B$ (I consider the real case for simplicity): $$B = \begin{bmatrix} P_{\img} & P_{\ker} \end{bmatrix}\begin{bmatrix} D & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} Q_{\img}^T \\ Q_{\ker}^T \end{bmatrix}.$$

Then set $x = Q_{\img}z_{\img} + Q_{\ker} z_{\ker}$ and plug this into the constraint $Bx=b$ you get $P_{\img}D z_{\img} = b$ and thus $z_{\img} = D^{-1}P_{\img}^Tb$. Note again that if $Bx=b$ has no solution the above would instead give you the solution to the projection $Bx = b^+ = P_{\img}P_{\img}^Tb$. Again denote $x_B = B^{+}_{\img}b = Q_{\img}D^{-1}P^T_{\img} b$, then $x = x_B + Q_{\ker}z_{\ker}$ and we can plug this into the minimisation: $$\min_{z_{\ker}} \frac{1}{2}(x_B+Q_{\ker}z_{\ker})^TL(x_B+Q_{\ker}z_{\ker}) - (x_B+Q_{\ker}z_{\ker})^Ta.$$ Taking the derivative w.r.t. $z_{\ker}$ and setting it equal to zero yields: $$Q_{\ker}^T\frac{L+L^T}{2}Q_{\ker}z_{\ker} - Q_{\ker}^T\left(a-\frac{L+L^T}{2}x_B\right) = 0.$$ Term $A = \frac{L+L^T}{2}$ (if $L$ is symmetric then $L=A$ like in your case). Then you have the compressed system: $$Q_{\ker}^TAQ_{\ker}z_{\ker} = Q_{\ker}^T\left(a-Ax_B\right).$$ here you have essentially projected your matrix $A$ onto the kernel of $B$. If a basis for the kernel of $B$ is tedious or computationally intensive to build, you could also use $Pr_{\ker(B)} = Q_{\ker}Q_{\ker}^* = I-Q_{\img}Q_{\img}^T = I-B^+B$ as a projection onto the kernel. Then you use the minimisation problem: $$\min_{y} \frac{1}{2}x_N^TPr_{\ker(B)}LPr_{\ker(B)}x_N - x_N^TPr_{\ker(B)}a,$$ which results in: $$Pr_{\ker(B)}APr_{\ker(B)}x_N = Pr_{\ker(B)}a.$$ To satisfy $Bx=b$ you also need to add $x_B = B^{+}b$ to this in order to get the solution $x = x_N+x_B$. It is of course clear that $x$ satisfies the compressed equation since $Pr_{\ker(B)}x_B = Pr_{\ker(B)}B^+b = 0$ because $Pr_{\ker(B)}B^+=0$ ($B^+b$ gives you a vector in the orthogonal space to $\ker(B)$).

If the projection is not cheap to compute, however, it may still be more efficient to just solve: $$\begin{bmatrix} A & B^T \\ B & 0 \end{bmatrix} \begin{bmatrix} x \\ \ell \end{bmatrix} = \begin{bmatrix} a \\ b \end{bmatrix},$$ for example with SYMMLQ, the biconjugate gradient stabilized solver, the conjugate residual solver, GMRES, or some other solver that works on symmetric indefinite systems. Note also that the minimisation formulation has a solution (unless $Bx=b$ has an empty solution set) even when the constrained system formulation doesn't have a solution.

System General Constrained Case

Here's a setting where I treat $A$ and $B$ as if they are "equally important".

$$\begin{bmatrix} A \\ B \end{bmatrix} x = \begin{bmatrix} a \\ b \end{bmatrix}.$$ Note that the solution of this problem would just be the intersection of two affine hyperplanes.

Suppose you have the singular value decomposition of $A$ and $B$: $$A = \begin{bmatrix} U_{\img} & U_{\ker} \end{bmatrix}\begin{bmatrix} \Sigma & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} V_{\img}^* \\ V_{\ker}^* \end{bmatrix}, \quad B = \begin{bmatrix} P_{\img} & P_{\ker} \end{bmatrix}\begin{bmatrix} D & 0 \\ 0 & 0 \end{bmatrix} \begin{bmatrix} Q_{\img}^* \\ Q_{\ker}^* \end{bmatrix}.$$

Then you can write $x = V_{\img}y_{\img} + V_{\ker} y_{\ker}$ and use that to solve $Ax=a$: $$Ax = A(V_{\img}y_{\img} + V_{\ker} y_{\ker}) = U_{\img}\Sigma y_{\img} = a \implies y_{\img} = \Sigma^{-1}U_{\img}^*a.$$ Be aware that if $a$ is not in the range of $A$ (i.e. if $Ax = a$ has no solution), then the above would actually give you the solution w.r.t. the projection of $a$ onto the range of $A$: $a_A = U_{\img}U_{\img}^*a$.

We can do the same thing for $B$ by setting $x = Q_{\img} z_{\img} + Q_{\ker} z_{\ker}$: $$Bx = P_{\img}D z_{\img} = b \implies z_{\img} = D^{-1}P_{\img}^* b.$$ A similar thing applies that the above would give you the solution wrt the projection $b_B = P_{\img}P_{\img}^* b$ if $Bx=b$ has no solution.

In either case the above equations did not set constraints on $y_{\ker}$ and $z_{\ker}$. We can now write equations for those also: \begin{align} V_{\img} y_{\img} + V_{\ker} y_{\ker} &= x = Q_{\img} z_{\img} + Q_{\ker} z_{\ker}\\ \begin{bmatrix} V_{\ker} & Q_{\ker} \end{bmatrix} \begin{bmatrix} y_{\ker} \\ -z_{\ker} \end{bmatrix} &= \underbrace{-V_{\img}\Sigma^{-1}U_{\img}^*a}_{-x_A}+ \underbrace{Q_{\img}D^{-1}P_{\img}^*b}_{x_B} = -x_A+x_B \end{align}

From the above it follows that if $\ker(A)\perp \ker(B)$ then a solution exists only if $x_A \in \ker(B)$ and $x_B\in \ker(A)$. If such a solution exists then it is unique and can be computed as follows: \begin{align} \begin{bmatrix} y_{\ker} \\ -z_{\ker} \end{bmatrix} &= \left(\begin{bmatrix} V_{\ker} & Q_{\ker} \end{bmatrix}^*\begin{bmatrix} V_{\ker} & Q_{\ker} \end{bmatrix}\right)^{-1}\begin{bmatrix} V_{\ker} & Q_{\ker} \end{bmatrix}^* (x_B-x_A) \\ &= \begin{bmatrix} I & V_{\ker}^*Q_{\ker} \\ Q_{\ker}^*V_{\ker} & I \end{bmatrix}^{-1} \begin{bmatrix} V_{\ker} & Q_{\ker} \end{bmatrix}^* (x_B-x_A) \\ &\stackrel{\ker(A)\perp\ker(B)}{=} \begin{bmatrix} I & 0 \\ 0 & I \end{bmatrix}^{-1} \begin{bmatrix} V_{\ker} & Q_{\ker} \end{bmatrix}^* (x_B-x_A) \\ &=\begin{bmatrix} V_{\ker}^*(x_B-x_A) \\ Q_{\ker}^*(x_B-x_A) \end{bmatrix} = \begin{bmatrix} V_{\ker}^*x_B \\ -Q_{\ker}^*x_A \end{bmatrix} \\ &= \begin{bmatrix} V_{\ker}^*Q_{\img}D^{-1}P_{\img}^*b \\ -Q_{\ker}^*V_{\img}\Sigma^{-1}U_{\img}^*a \end{bmatrix}. \end{align} If we introduce the notation $A^{+}_{\img} = V_{\img}\Sigma^{-1}U_{\img}^*$, $B^{+}_{\img} = Q_{\img}D^{-1}P_{\img}^*$ then: \begin{alignat}{2} y_{\img} &= V_{\img}^*x_A = V_{\img}^*A^{+}_{\img}a, \quad &y_{\ker} &= V_{\ker}^*x_B = V_{\ker}^*B^{+}_{\img} b, \\ z_{\img} &= Q_{\img}^*x_B = Q_{\img}^*B^{+}_{\img}b,\quad &z_{\ker} &= Q_{\ker}^*x_A = Q_{\ker}^*A^{+}_{\img} a. \end{alignat} From the above it follows that $x = x_A + x_B$. If $\ker(A)\perp \ker(B)$ but $x_A \not\in \ker(B)$ or $x_B \not\in \ker(A)$, then $x = x_A+x_B = A^+a + B^+b$ is not a solution. It is equivalent to having taken the Moore-Penrose pseudoinverse solutions of the two problems and adding them up.

For the case where $\ker(A)\not\perp\ker(B)$ there are infinitely many solutions if $(x_B-x_A) \in \ker(A)+\ker(B)$. That means that insufficiently many constraints were supplied. It is clear that $y_{\ker} = V_{\ker}^*x_B$ and $z_{\ker} = Q_{\ker}^*x_A$ is a solution, but I don't think it is necessarily the Moore-Penrose pseudoinverse one.

  • 1
    $\begingroup$ Thank you very much for the detailed answer. I will read it carefully to understand the ideas. $\endgroup$ Feb 23 at 21:24
  • 1
    $\begingroup$ @BeniBogosel I had a typo which I fixed, namely $Pr_{\ker(B)} = I-Q_{\img}Q_{\img}^* = I-B^+B$. I had previously written $I-BB^+$ which is obviously wrong since it gives $I-P_{\img}P_{\img}^*$, sorry about that. This was in the general setting though, so I hope that it didn't cause any issue. $\endgroup$
    – lightxbulb
    Feb 26 at 13:58
  • $\begingroup$ It turns out that in the case I investigate, the simple conjugate gradient with initialization $x=0$ is faster than solutions that I tried previously. Thanks for the suggestion. $\endgroup$ Mar 25 at 16:37
  • $\begingroup$ @BeniBogosel Conjugate gradients is often times very efficient (especially if you precondition it properly). Even if the constraints were not from the kernel you can still use an equality constrained conjugate gradient method, especially if the projection on the nullspace of the constraints is efficient (e.g. either analytically known, or because you can construct an orthonormal basis cheaply, or using CGNR is cheap). See for example the paper: "On the Solution of Equality Constrained Quadratic Programming Problems Arising in Optimization", epubs.siam.org/doi/10.1137/S1064827598345667 $\endgroup$
    – lightxbulb
    Mar 25 at 16:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.