Solving linear system using Sherman-Morrison formula for 1000000x1000000 (7450.6GB) matrix

Let . Let A be the lower triangular matrix having 1's on and below the main diagonal.

We want to solve the following linear system:

by the Sherman-Morrison formula:

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

Suppose and suppose be the solution of be the solution of Then the solution of is given by

My question is how do I compute I know A is a unit lower triangular matrix so the formula is for and for and because all the entries below the main diagonal are 1 , for I know this is Forward substitution but how do I incorporate this in MATLAB?

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 an algorithm in computing if it is impossible to store A.

[Christine Tobler]

This sounds a bit like a homework problem. You should take a look at what inv(triu(ones(100))) looks like - you probably should never try to allocate that large of a matrix. Avoiding that seems to be the whole point of your problem here.

Answer 1

[Steven Lord]

Unless you have a machine with multiple terabytes of memory, you can't solve this directly.

Instead, use one of the iterative solvers. Most if not all of them have the capability to allow you to specify either an explicitly created matrix or a function handle that given a vector x as input computes A*x and/or A'*x.

You want to solve the system (A+u*v')*x = b. If we expand the left side we get (A*x)+(u*v'*x) = b. That first term, given the structure of A, is easy to compute without explicitly creating A (cumsum). I'll leave it to you to think about how to compute (u*v'*x) without explicitly creating the matrix u*v'.

[Steven Lord]

Actually, on second thought my first sentence there isn't quite correct. If you had a cluster of machines that you could use with MATLAB Parallel Server and Parallel Computing Toolbox you could use some of the capabilities for working with Big Data. The backslash operator \ is supported for use with distributed arrays.

But creating that large distributed array would still take time, so I'd still take a look at the iterative solvers first.