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

10 次查看（过去 30 天）
Purnjay Peshawaria
Purnjay Peshawaria 2020-6-30
评论： Christine Tobler 2020-6-30
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 个评论
Christine Tobler
Christine Tobler 2020-6-30
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Steven Lord
Steven Lord 2020-6-30
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 个评论
Steven Lord
Steven Lord 2020-6-30
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Operating on Diagonal Matrices 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by