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

Question

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}.$$

We are asked to compute:

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?

Answer 1

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

Answer 2

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$