# Solution of the linear system using Sherman-Morrison formula for 1000000x1000000 (7450.6GB) matrix using MATLAB

Let $$n = 10^6.$$ Let $$A \in \mathbb{R}^{n\times n}$$ be the lower triangular matrix having 1's on and below the main diagonal.

We want to solve the following linear system: $$(A + uv^T)x = b$$

by the Sherman-Morrison formula: $$(A+uv^T)^{-1} = A^{-1}-\frac{A^{-1}uv^TA^{-1}}{1+v^TA^{-1}u}.$$

1. $$A^{-1}b$$
2. $$1+ v^TA^{-1}u$$
3. $$v^TA^{-1}b$$
4. Solution of the linear system using Sherman-Morrison formula

My attempt:

% initialse n
n = 1e6;
% generate random vectors u,v,b
rng(1);
u = randn(n,1);
v = randn(n,1);
b = randn(n,1);
% create lower triangular matrix having 1's on and below the main diagonal
A = tril(ones(n,n));


I'm getting the following error:

Error using ones Requested 1000000x1000000 (7450.6GB) array exceeds maximum array size preference. Creation of arrays greater than this limit may take a long time and cause MATLAB to become unresponsive. See array size limit or preference panel for more information.

I need help with storing A and solving the linear system.

I have spent some time reading about Sherman-Morrison formula. Here is what I have understood:

Suppose $$det(A) \neq 0$$ and $$det(A + uv^T) \neq 0$$ and suppose $$\mathbf x = \mathbf x^* \in \mathbb {R}^n$$ be the solution of $$A\mathbf x = \mathbf b, \mathbf y = \mathbf y^* \in \mathbb {R}^n$$ be the solution of $$\mathbf A \mathbf y = \mathbf u.$$ Then the solution of $$(A + \mathbf u \mathbf v^T)\mathbf x = \mathbf b$$ is given by $$\mathbf x = \mathbf x^* - \frac{\mathbf v^T \mathbf x^*}{1+\mathbf v^T \mathbf y^*} \mathbf y^*.$$

But again my question is how do I compute $$A^{-1} \mathbf b.$$ I know A is a unit lower triangular matrix so the formula is $$a_{ij} = 0$$ for $$1 \leq i < j \leq n$$ and $$a_{ii} = 1$$ for $$1 \leq i \leq n$$ and because all the entries below the main diagonal are 1 , $$a_{ij} = 1$$ for $$1 \leq j < i \leq n.$$ I know this is Forward substitution but how do I incorporate this is MATLAB?

• This looks a lot like it is an exercise given to you especially so that you'd stumble onto this problem and need to find a solution. Are you sure it is a good idea to deprive yourself of the learning experience? Jun 27 '20 at 6:20

As Federico has mentioned, you probably don't want to deprive yourself of the learning experience. I'll just give you a small nudge in the right direction.

You will never be able to store $$A$$. You also won't be able to store $$(A+uv^T)^{-1}$$. However, you don't really need to. You can easily write down a formula for each of the entries in $$A$$.

Instead of relying on matlab to solve the system for you, look up some algorithms for solving triangular matrices. Everytime you need to use $$A_{i,j}$$ in some formula, substitute it by either $$1$$ or $$0$$.

You can apply similar methods to solve the other questions.

edit, here is some code illustrating forward substitution, you should be able to extrapolate from there:

n = 10;
A = tril(ones(n,n));
b = rand(n,1);
x = zeros(n,1);

for k1=1:n
x(k1) = b(k1);
for k2 = 1:k1-1
x(k1) = x(k1) - A(k1,k2)*x(k2);
end
x(k1) = x(k1)/A(k1,k1);
end

• I have modified my question. Please have a look. I know A is a unit lower triangular matrix and we can solve it by Forward substitution. But how do I implement this in MATLAB?
– user36184
Jun 30 '20 at 6:42
• i've edited my answer. BTW, don't you have a TA that you can contact for this sort of thing? Jun 30 '20 at 9:20

The inverse of $$D$$ will form a matrix such that

The diagonal of the inverse of $$D$$ is $$1$$ (for $$i=1,2,...n$$) $$D(i,i) = 1$$

and for $$i=1,2,...n-1$$ $$D(i,i+1) = -1$$

and the rest of the values of the inverse of $$D$$ are zero.

Example (using $$n=5$$)

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

So now you can use the values of $$D$$ and $$D^{-1}$$ to compute the terms involved in the equations and solve for x without storing $$D$$