7

2

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)?

flag
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 at 12:27
3 
 If you don't need the exact solution but some approximation, you can use an iterative method like Gauss-Seidel. – Brendan McKay Dec 7 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 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 at 16:03

2 Answers

6

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.

link|flag
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 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 at 12:21
3

Being able to get elements of the matrix isn't very useful (particularly if you don't know where the nonzero elements of the matrix are without checking.)

Iterative methods can be useful if you have the ability compute matrix vector products $Mx$. You haven't said whether this is possible.

It seems quite likely that there's some special structure to your particular problem that would make it possible to simplify this computation. You haven't told us anything about where the system of equations comes from- perhaps if you explained this in some detail we could suggest ways to proceed.

link|flag
3 
 Being able to get elements of the matrix is useful to compute matrix vector products, for instance. Am I missing something? – Emil Jeřábek Dec 7 at 14:48
As already said above, being able to get matrix elements is enough to compute matrix vector products. The whole matrix consists mostly of non-zero elements, so it's not useful to search for any zero ones. As for structure, the matrix is an effective Hamiltonian, so what I can say for sure is only that it's Hermitian. As for origin of the system of equations, I'm trying to use inverse iteration with shift to find specific eigenvectors of this matrix, so Y is current approximation of eigenvector, and X is next step approximation. – 10110111 Dec 7 at 15:54
Let me clarify what I meant here- "being able to get an arbitrary element M(i,j) at little cost" isn't very useful. If you don't know where the nonzero elements are in the matrix, then you have to check every single one to find the nonzeros. If you do happen to know where the nonzero elements are, and you can compute them quickly, then you could use this as a way to do matrix vector multiplications in an iterative method. – Brian Borchers Dec 16 at 14:22
If the matrix isn't sparse, and the cost of getting individual matrix entries is large compared to the cost of accessing an element of a matrix stored in conventional dense matrix form, then iterative methods are going to be horribly slow in practice. – Brian Borchers Dec 16 at 14:25

Your Answer

Get an OpenID
or

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