# Find the inverse of a submatrix of a given matrix

I have a $A$ matrix $4 \times 4$. By delete the first row and first column of $A$, we have a matrix $B$ with sizes $3 \times 3$. Assume that I have the result of invertible A that denote $A^{-1}$ before. I want to use the result to calculate the inverse of $B.$ Is it possible? Have any method to do it? Thanks

Let consider a simple case that A is invertible then B also invertible.

• There's no reason for $B$ to be invertible. Commented Apr 23, 2015 at 12:37
• Let assume that B is invertible
– John
Commented Apr 23, 2015 at 13:06

Note that $A$ differs from $\begin{bmatrix} 1 & 0 \\ 0 & B \end{bmatrix}$ by at most a rank two "update". Therefore formulas which "update" the matrix inverse $A^{-1}$ can be used to give an inverse of $B$:

$$\begin{bmatrix} 1 & 0 \\ 0 & B \end{bmatrix}^{-1} = \begin{bmatrix} 1 & 0 \\ 0 & B^{-1} \end{bmatrix}$$

A rank one update to $A$ corresponds to a rank one update to $A^{-1}$ by the Sherman-Morrison formula:

$$(A + uv^T)^{-1} = A^{-1} - \frac{A^{-1}uv^T A^{-1} }{1 + v^T A^{-1} u}$$

The more general Woodbury formula could be used here to give the inverse of the rank two update to $A$ directly as a rank two update to $A^{-1}$.

$$\left(A+UCV \right)^{-1} = A^{-1} - A^{-1}U \left(C^{-1}+VA^{-1}U \right)^{-1} VA^{-1}$$

The parameter $C$ above is redundant in the sense that we can take it to be the identity (and this is often done) choosing $U,V$ to "absorb" any factor $C$. In the case at hand:

$$A = \begin{bmatrix} a & v_1 & v_2 & v_3 \\ u_1 & & & \\ u_2 & & B & \\ u_3 & & & \end{bmatrix}$$

we can take $U = \begin{bmatrix} 1 & 0 \\ 0 & u_1 \\ 0 & u_2 \\ 0 & u_3 \end{bmatrix}$ and $V = \begin{bmatrix} 0 & v_1 & v_2 & v_3 \\ 1&0&0&0 \end{bmatrix}$, so that:

$$UV = \begin{bmatrix} 0 & v_1 & v_2 & v_3 \\ u_1 & 0 & 0 & 0 \\ u_2 & 0 & 0 & 0 \\ u_3 & 0 & 0 & 0 \end{bmatrix}$$

Then $A - UV = \begin{bmatrix} a & 0 \\ 0 & B \end{bmatrix}$ and provided $a \neq 0$ the update can be done by inverting $(I - VA^{-1}U)$, a $2\times 2$ matrix.

For the small sizes of matrices $A$ and $B$ (resp. $4\times 4$ and $3\times 3$) such computations may promise little improvement over a direct computation of $B^{-1}$. As the Comments indicated interest in removing multiple $k$ columns and rows from larger $n\times n$ matrices, the Woodbury formula can be used for this. However one needs to keep in mind the cost of inverting the $2k\times 2k$ matrix $I + VA^{-1}U$ when $2k$ is a significant fraction of original size $n$.

• It is very good answer. Actually, I am working in very large matrix such as $1000 \times 1000$ or more smaller is $500 \times 500$. I also saw a paper used Sherman-Morrion to resolve my problem. Hence, when I see your answer again. It make to me more information. Could you know any more faster or improvement of Sherman-Morrison method? Thanks for your help
– John
Commented Apr 24, 2015 at 17:48
• So, $A$ is perhaps $1000\times 1000$ and $B$ is $500\times 500$? The update approach is attractive when the submatrix $B$ is nearly the same size as $A$, but not in cases where the sizes are so different. The first thing to ask is whether a matrix inverse is really needed. In most numerical methods, being able to solve linear systems is more expedient than computing inverses. Commented Apr 24, 2015 at 17:56
• No no no. I just give two example of matrix $A$ size. One for large matrix $1000 \times 1000$. And one is smaller that $500 \times 500$. The matrix $B$ is created from $A$ and deleted some cols and rows. Number of rows and cols around 10% number of cols or rows in $A$
– John
Commented Apr 24, 2015 at 18:00
• And I known one common method in my area is belief propagation. However, it only well apply for sparse matrix. But my matrix is not it.
– John
Commented Apr 24, 2015 at 18:07
• I've added some material about doing all the rows and columns in a single update. However my instinct is that doing the rows OR the columns ought to suffice, if what you are after is merely $B^{-1}$. I'll have to find time to explain this thought. Commented Apr 28, 2015 at 15:40

http://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion

or

http://en.wikipedia.org/wiki/Block_matrix_pseudoinverse

The pseudoinverse will have to be used in case B is not invertible. For instance the case in the other answer by 5xum.

• In particular you may want to familiarize yourself with the "Schur complement" en.wikipedia.org/wiki/Schur_complement which appears in many contexts. Commented Apr 23, 2015 at 13:14
• Thank mathreadler. It is very interesting answer. However, is it possible to apply Schur for inverst matrix
– John
Commented Apr 23, 2015 at 13:37
• The blockwise inversion works only if all submatrices are blocks. Commented Dec 23, 2019 at 20:35

Block-inverse and Schur-complement were both mentioned but not developed into a closed form expression, so perhaps the following is still useful.

Defining a unit vector $$r$$ that selects the single row and column to be deleted from the matrix, and the $$n \times n-1$$ complement matrix $$R$$ such that $$R R^T + r r^T = I$$, we can block-decompose a matrix as $$M = R A R^T + R B r^T + r C R^T + r D r^T,$$ where $$A = R^T M R$$, $$B = R^T M r$$, $$C = r^T M R$$ and $$D = r^T M r$$. If $$D$$ is nonzero and the Schur-complement $$M/D$$ is invertible we can similarly decompose the inverse as $$M^{-1} = R (M/D)^{-1} R^T + R \dots r^T + r \dots R^T + r \dots r^T.$$

Isolating the first block (using $$r^T R = 0$$ and $$R^T R = I$$) we obtain $$(R^T M^{-1} R)^{-1} = M / D = A - B D^{-1} C = R^T \left( M - \frac {M r r^T M} {r^T M r} \right) R.$$

Substituting $$M = P^{-1}$$ we obtain an expression for the inverse of a submatrix of $$P$$ in terms of $$P^{-1}$$: $$(R^T P R)^{-1} = R^T \left( P^{-1} - \frac {P^{-1} r r^T P^{-1}} {r^T P^{-1} r} \right) R.$$

Here I discuss how to generalize hardmath's (excellent) answer as follows:

• $$k$$ rows and $$k$$ columns are removed instead of just one.
• The rows and columns may be at any location in the matrix (not necessarily the front, and not necessarily in a contiguous block)

I came across this question because I needed to solve this more general problem for a numerical method I am working on. Conceptually, all the key ideas are in hardmath's answer, but I found it still took a bit of work to do the generalization to a level of detail where I could implement it on a computer efficiently. I am posting my findings here to save future people the trouble. I also include python/numpy code which anyone may use as they see fit.

Let $$\mathcal{R}^g$$ be an index set of "good" rows of $$A$$, and let $$\mathcal{R}^b$$ be the complementary index set of "bad" rows of $$A$$, so that $$\mathcal{R}^g \cup \mathcal{R}^b = \{1,2,\dots,N\}$$, where $$A$$ is an $$N \times N$$ matrix.

Likewise, let $$\mathcal{C}^g$$ be an index set of "good" columns of $$A$$, and $$\mathcal{C}^b$$ be the complementary set of "bad" columns of $$B$$.

We seek to solve the linear system $$A\left(\mathcal{R}^g,\mathcal{C}^g\right) ~x = b, \qquad (1)$$ where $$A\left(\mathcal{R}^g,\mathcal{C}^g\right)$$ denotes the submatrix of $$A$$ consisting of the "good" rows and columns. (The same notation as how you select a submatrix in Matlab)

For example, if $$\mathcal{R}^g=(1,3)$$, $$\mathcal{C}^g=(5,4)$$, and $$N=5$$, then the matrix in equation (1) becomes $$A\left(\mathcal{R}^g,\mathcal{C}^g\right) = A\left((1,3),(5,4)\right) = \begin{bmatrix} A_{15} & A_{14} \\ A_{35} & A_{34} \end{bmatrix},$$ and the complementary "bad" index sets are $$\mathcal{R}^b=(2,4,5)$$ and $$\mathcal{C}^b=(1,2,3)$$. (the ordering of the bad index sets is not important so long as one keeps it consistent)

Now consider the matrix $$\widetilde{A}$$ which is the same as $$A$$, except matrix entries corresponding to interactions between "good" and "bad" index sets are set to zero. I.e., \begin{align} \widetilde{A}\left(\mathcal{R}^g,\mathcal{C}^g\right) &= A\left(\mathcal{R}^g,\mathcal{C}^g\right) \\ \widetilde{A}\left(\mathcal{R}^b,\mathcal{C}^b\right) &= A\left(\mathcal{R}^b,\mathcal{C}^b\right) \\ \widetilde{A}\left(\mathcal{R}^g,\mathcal{C}^b\right) &= 0 \\ \widetilde{A}\left(\mathcal{R}^b,\mathcal{C}^g\right) &= 0. \end{align} Also, let $$\widetilde{b}$$ be the vector with $$\widetilde{b}(\mathcal{R}^g)=b$$, and $$\widetilde{b}(\mathcal{R}^b)=0$$, where $$\widetilde{b}(\mathcal{R}^g)$$ denotes the entries of $$\widetilde{b}$$ corresponding to indices in $$\mathcal{R}^g$$, and likewise for $$\widetilde{b}(\mathcal{R}^b)$$.

Since these cross interactions in $$\widetilde{A}$$ are zero, if we solve $$\widetilde{A} \widetilde{x} = \widetilde{b}, \qquad (2)$$ then $$x\left(\mathcal{C}^g\right) = x.$$ In other words, the solution of (1) is found by extracting the $$\mathcal{C}^g$$ entries from the solution of (2).

Now, like in hardmath's answer, $$\widetilde{A}$$ can be written as a low rank update of $$A$$ in the form $$\widetilde{A} = A - UV,$$ and system (2) can be solved via the Woodbury formula: \begin{align} \widetilde{x} = \widetilde{A}^{-1}\widetilde{b} &= \left(A - UV\right)^{-1}\widetilde{b} \\ &= A^{-1}\widetilde{b} + A^{-1}U\left(I_{2k \times 2k} - VA^{-1}U\right)^{-1}VA^{-1}\widetilde{b}. \end{align} But here $$U$$ is a $$N \times 2k$$ matrix, $$V$$ is a $$2k \times N$$ matrix, and $$I_{2k \times 2k}$$ is the $$2k \times 2k$$ identity matrix, where $$k$$ is the number of "bad" rows/columns (size of $$\mathcal{R}^b$$ and $$\mathcal{C}^b$$).

The matrices $$U$$ and $$V$$ are given as follows: \begin{align} U\left(\mathcal{R}^g, (1,\dots,k)\right) &= 0 \\ U\left(\mathcal{R}^b, (1,\dots,k)\right) &= I_{k\times k} \\ U\left(\mathcal{R}^g, (k+1,\dots,2k)\right) &= A\left(\mathcal{R}^g, \mathcal{C}^b\right) \\ U\left(\mathcal{R}^b, (k+1,\dots,2k)\right) &= 0 \end{align} \begin{align} V\left((1,\dots,k), \mathcal{C}^g\right) &= 0 \\ V\left((1,\dots,k), \mathcal{C}^b\right) &= A\left(\mathcal{R}^b, \mathcal{C}^g\right) \\ V\left((k+1,\dots,2k), \mathcal{C}^g\right) &= I_{k\times k} \\ V\left((k+1,\dots,2k), \mathcal{C}^b\right) &= 0. \end{align}

Here is some python/numpy code demonstrating how to do this.

First, building the matrices $$U$$ and $$V$$ and verifying that they are correct:

import numpy as np

N = 102
A = np.random.randn(N, N)

bad_rows = [2,   1, 37, 15]
bad_cols = [80, 60, 40, 55]

A_tilde_true = A.copy()

U = np.zeros((A.shape[0], 2*k))

V = np.zeros((2*k, A.shape[1]))

A_tilde = A - np.dot(U, V)
err_UV = np.linalg.norm(A_tilde_true - A_tilde) / np.linalg.norm(A_tilde_true)
print('err_UV=', err_UV)


which outputs: "err_UV= 0.0"

And now using $$U$$ and $$V$$ with the Woodbury formula to solve a linear system with the submatrix:

b = np.random.randn(N-k)
x_true = np.linalg.solve(A[np.ix_(good_rows, good_cols)], b)

Z = np.linalg.solve(A, U) # Z = inv(A)*U
C = np.eye(2*k) - np.dot(V, Z) # C = I - V*inv(A)*U

b_tilde = np.zeros(N)
b_tilde[good_rows] = b

x_tilde = np.linalg.solve(A_tilde, b_tilde) # x0 = inv(A)*b
x_tilde += np.dot(Z, np.linalg.solve(C, np.dot(V, x_tilde))) # dx = inv(A)*U*inv(C)*V*x0
x = x_tilde[good_cols]

err_woodbury = np.linalg.norm(x - x_true) / np.linalg.norm(x_true)
print('err_woodbury=', err_woodbury)


which outputs "err_woodbury= 9.049341577063206e-14" (i.e., zero, up to floating point rounding errors).

You can find an example in which $B$ is not invertible:

$$A=\begin{bmatrix}0 & 1 & 0 & 0\\ 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 1 \end{bmatrix}$$

meaning that what you want to do is impossible.

• A^2 = I so A is it's own inverse. Maybe you mean that it is not diagonalizable over it's scalar field of elements. Commented Apr 23, 2015 at 12:51
• Oh, disregard what I wrote. I misread what you wrote. You are correct that B is not invertible, I read it as if you meant A. Commented Apr 23, 2015 at 12:53
• @5xum: let assume that B is invertible. Just simple case
– John
Commented Apr 23, 2015 at 13:07