I'm a beginner in linear algebra, so my notation and fundamentals might be sloppy.

Suppose that the basis for the null space of a matrix $\mathbf{A}$, size $N \times N$, is a vector of length $N$ with all entries equal to $1$, i.e. \begin{align} \mathrm{ker}(\mathbf{A}) = \mathrm{span}(\mathbf{1}) \end{align} where $\mathbf{1}$ is the aforementioned column vector of ones, with size $N \times 1$. Therefore, \begin{align} \mathbf{A}\mathbf{1} = \mathbf{0}. \end{align} For context, $\mathbf{A}$ represents a boundary element method discretization of the Laplace equation, for the interior Neumann problem for electrostatics. The fact that $\mathbf{1}$ is a basis for the null space represents the physical fact that the electrostatic potential can only be determined up to a constant for the interior Neumann problem. The actual value of the potential must then be determined using the exterior problem.

I'm interested in solving the system \begin{align} \mathbf{A}\mathbf{x} = \mathbf{b},\tag{1} \end{align} eventually with an iterative solver. I can write \begin{align} \mathbf{x} = \mathbf{1}x_0 + \mathbf{P}\widetilde{\mathbf{x}},\tag{2} \end{align} where $\mathbf{P}$ is some matrix which, I believe, should be a basis for $\mathrm{range}(\mathbf{A})$. Knowing $\mathbf{P}$ may help solve the system in the following way: the idea is to split the unknowns into a constant offset $x_0$ and the remaining non-constant part $\widetilde{\mathbf{x}}$. Then, one can solve the interior problem to determine $\widetilde{\mathbf{x}}$. Finally, by coupling (1) to the exterior problem, one should be able to determine $x_0$ separately. Finally, we can use (2) to obtain the final solution $\mathbf{x}$.

My two questions are: (1) Is this a valid approach? (2) If so, is there an easy way to determine $\mathbf{P}$ (or even just the product of $\mathbf{P}$ with a vector) given the above information, which does not involve converting $\mathbf{A}$ to echelon form?

Note: I also found literature on the use of the "deflation" approach, in which the system of equations is written as \begin{align} \mathbf{A}\left(\mathbf{I} - \mathbf{1}(\mathbf{1}^T\mathbf{1})^{-1}\mathbf{1}^T\right)\mathbf{x} = \mathbf{b}.\tag{3} \end{align} My understanding is that the solution space is projected to the span of all vectors which are orthogonal to the null space, or in simple (overly simple?) terms, the null space is "avoided / cancelled". I feel like this method should be related to my approach of trying to find $\mathbf{P}$. I realize that the latter approach might be better, but I still want to pursue the former approach of trying to determine $\mathbf{P}$, for the sake of understanding and reconciling these two approaches, which to me seem somehow equivalent.

Edit: another reason I want to make the former method work is that in the latter method, the matrix $\mathbf{A}\left(\mathbf{I} - \mathbf{1}(\mathbf{1}^T\mathbf{1})^{-1}\mathbf{1}^T\right)$ is singular - it is rank deficient with one very tiny singular value. This worries me, because I want to eventually solve the system with an iterative solver. However, I am open to ideas for dealing with this singularity?


1 Answer 1


I think I've figured out the answer to my first question, in case it is helpful to anyone: the former approach is valid, because it is basically equivalent to the latter approach.

Consider again the system \begin{align} \mathbf{A}\mathbf{x} = \mathbf{b}\tag{1}\label{eq1} \end{align} In the latter deflation approach, we instead solve the system \begin{align} \mathbf{A}\left(\mathbf{I} - \mathbf{1}(\mathbf{1}^T\mathbf{1})^{-1}\mathbf{1}^T \right)\mathbf{x}_r = \mathbf{b}\tag{2}\label{eq2} \end{align} for $\mathbf{x}_r$. According to the former approach of splitting the unknowns, we instead want to write \begin{align} \mathbf{A}\left(\mathbf{1}x_0 + \mathbf{P}\widetilde{\mathbf{x}}\right) = \mathbf{b}\tag{3}\label{eq3} \end{align} Comparing (\ref{eq2}) and (\ref{eq3}), \begin{align} \mathbf{1}x_0 + \mathbf{P}\widetilde{\mathbf{x}} &= \left(\mathbf{I} - \mathbf{1}(\mathbf{1}^T\mathbf{1})^{-1}\mathbf{1}^T \right)\mathbf{x}_r\\ &= \mathbf{x}_r + \mathbf{1}\underbrace{-(\mathbf{1}^T\mathbf{1})^{-1}}_{\text{constant}}\underbrace{\mathbf{1}^T\mathbf{x}_r}_{\text{constant}}\\ &= \mathbf{x}_r + \mathbf{1}x_0\tag{4}\label{eq4} \end{align} (The last step is somewhat by intuition / hypothesis). Rearranging (\ref{eq4}), \begin{align} \mathbf{1}x_0 + \mathbf{P}\widetilde{\mathbf{x}} &= \mathbf{x}_r + \mathbf{1}x_0\\ \implies \mathbf{x}_r &= \mathbf{P}\widetilde{\mathbf{x}}\tag{5}\label{eq5} \end{align} Now notice from (\ref{eq3}) that \begin{align} \mathbf{A}\left(\mathbf{1}x_0 + \mathbf{P}\widetilde{\mathbf{x}}\right) &= \mathbf{b}\\ \implies \underbrace{\mathbf{A}\mathbf{1}x_0}_{0 \text{ (null space)}} + \mathbf{A}\mathbf{P}\widetilde{\mathbf{x}}&= \mathbf{b}\\ \implies \mathbf{A}\mathbf{P}\widetilde{\mathbf{x}}&= \mathbf{b}\tag{6}\label{eq6} \end{align} Therefore, solving (\ref{eq2}) for $\mathbf{x}_r$ is equivalent to solving (\ref{eq6}) for the vector $\mathbf{P}\widetilde{\mathbf{x}}$; in both cases, the desired vector of unknown $\mathbf{x}$ can be obtained as a simple post-processing step.

I still don't know the answer to my second question.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .