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

32 visualizaciones (últimos 30 días)
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.
  1 comentario
Christine Tobler
Christine Tobler el 30 de Jun. de 2020
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.

Iniciar sesión para comentar.

Respuestas (1)

Steven Lord
Steven Lord el 30 de Jun. de 2020
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'.
  1 comentario
Steven Lord
Steven Lord el 30 de Jun. de 2020
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.

Iniciar sesión para comentar.

Categorías

Más información sobre Arithmetic Operations en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by