2
$\begingroup$

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?

$\endgroup$
  • 12
    $\begingroup$ This looks a lot like it is an exercise given to you especially so that you'd stumble onto this problem and need to find a solution. Are you sure it is a good idea to deprive yourself of the learning experience? $\endgroup$ – Federico Poloni Jun 27 at 6:20
6
$\begingroup$

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
| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ I have modified my question. Please have a look. I know A is a unit lower triangular matrix and we can solve it by Forward substitution. But how do I implement this in MATLAB? $\endgroup$ – Rhombus Jun 30 at 6:42
  • 2
    $\begingroup$ i've edited my answer. BTW, don't you have a TA that you can contact for this sort of thing? $\endgroup$ – Thijs Steel Jun 30 at 9:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.