# How to solve a system of linear equations without storing the matrix?

I have a procedurally defined Hermitian matrix $M$, i.e. I can get any matrix element by calling a black box function (e.g. a library function), and a vector $Y$. And I have to solve a system of linear equations:
$M\cdot X=Y$.
But $Y$ is such that having $n$ elements, it takes about half of available RAM, another half would be for $X$, so if I try to store $M$, it'll take $n^2$ space, i.e. even if I double the RAM space, this will still place $n-2$ matrix columns/rows into swap.
At the same time, all the algorithms which I've found need saving large amounts ($\geq n^2$) of data while solving the system.

So, the question: are there any algorithms to solve systems of linear equations which don't require me to store the matrix, and still are fast enough (maybe not $O(n^3)$, but at least not much slower)?

-
In principle, determinants, and therefore solutions of non-singular linear systems, are computable in $\mathrm{NC}^2$, and therefore in space $O(\log^2n)$ (not including the input and output). However, I rather doubt such algorithms would be practical. In particular, they need time $n^{O(\log n)}$, which most likely wouldn’t count as “not much slower”. –  Emil Jeřábek Dec 7 '12 at 12:27
If you don't need the exact solution but some approximation, you can use an iterative method like Gauss-Seidel. –  Brendan McKay Dec 7 '12 at 12:40
I assume reading parts of the matrix, modify them, and write back is out of the question, since this is essentially as having a huge swap...? –  Per Alexandersson Dec 7 '12 at 15:16
@Per Alexandersson Of course, this is not an option. Even if HDD/SSD speed were comparable with RAM speed, it'd take petabytes of space for about a gig of RAM. –  10110111 Dec 7 '12 at 16:03

This looks like a situation where the Kaczmarz method could work.

What you do to maintain an approximate solution and then project cyclically onto the hyperplanes which are given by the $k$-th equation. More precisely: If you have the $m$-the iterate $X^m$ and use the $k$-th equation, then you have the next iterate $$X^{m+1} = X^M + \frac{Y_k - a_k^T\cdot X^m}{\|a_k\|^2}a_i$$ where $a_k$ is the $k$-th row of $A$ and $Y_k$ is the $k$-the entry of $Y$.

Hence, you only need one row of $A$, one entry of $Y$ and the current iterate $X^m$ to perform one iteration, i.e. $2n+1$ space. Also the iteration complexity is very low (it's $\mathcal{O}(n)$) but you usually need a lot of iterations.

This methods is widely used in discrete tomography and is also an instance of the "Projection onto convex sets" method.

Recently it has been shown by Strohmer and Vershynin that a randomized version of this method has favorable convergence properties (when you pick each column with a probability proportional to its norm). Also "block iterative" versions work, i.e. you take a hyperplane of higher codimension to project on. So, if you have some memory left, you could also take some more rows of $A$ at once... See, e.g. here.

-
Well, indeed converges rather slowly. Though, I've not tried randomized version yet.<br /> Do I understand correctly that I can use this method without requirements of matrix being real or diagonally dominant as is for Gauss-Seidel? –  10110111 Dec 9 '12 at 11:23
Granted, convergence can be slow (in terms of iteration count and computational effort - its advantage is low memory). You don't need any requirements for the matrix (for the complex case adjust the projection accordingly). In fact you could also apply the method to rectangular systems. It converges to some solution for the underdetermined case (and the minimum-norm solution if initialized with zero). In the overdetermined case you need to stop at some point as you'll see that the residual $\|AX-Y\|$ is not decreasing anymore. –  Dirk Dec 9 '12 at 12:21
Iterative methods can be useful if you have the ability compute matrix vector products $Mx$. You haven't said whether this is possible.