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:
- $A^{-1}b$
- $1+ v^TA^{-1}u$
- $v^TA^{-1}b$
- 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?