Iterative solver for matrix division

I have two huge matrices $A$ and $B$. I am trying to find some iterative solvers like bcg or lsqr in Matlab.

I mean if my matrix $A$ is sparse and $n\times n$ and $B$ is a column vector of size $n\times 1$, I can use lsqr and bcg to solve the equation $Ax=B$.

Now suppose I need to solve $XD=C$, so I need to calculate $CD^{-1}$ where both $C$ and $D$ are huge matrices. If I use matlab's C/D operation directly, it consumes lots of memory and crashes. Are there any iterative solvers for this operation instead of using the forward slash operator directly?

-
How big is a "huge" matrix? – Chris Taylor Oct 11 '12 at 11:44
Around 4 million by 4 million – user34790 Oct 11 '12 at 11:50
If you can transpose easily, solve $D^TX^T=C^T$ by standard techniques. – Daryl Oct 11 '12 at 11:56
But methods like lsqr, bcg require the matrix B in the equation Ax=B to be a column matrix of size nx1. In my case the B matrix will not be a column matrix.So what should I do – user34790 Oct 11 '12 at 11:58
Have you tried writing a loop which solves each of the columns of B individually? For example result = sparse(N,N); for i = 1:N; result(:,i) = A\B(:,i); end. – Chris Taylor Oct 11 '12 at 12:06

The function GMRES offers the best speed, though I think QMR uses less memory. Otherwise the lu function allows you to recompose the matrix into an upper and a lower matrix like so: [L,U,P] = lu(A); x = U(L(P*b)); Where A*x = b

-
LU isn't that computationally feasible for huge problems... – Daryl Oct 29 '12 at 23:08
Right, but I think it intersects CG at around N=8000, so it's alright for small problems. Big ones I'd go with QMR or GMRES. – Rorrik Oct 31 '12 at 18:28

Let's frame your problem as a Sylvester equation:

$$AX+XD - C = 0$$

with $A = 0$. We can solve Sylvester equations using lyap.

As per my comment, the following will work:

n = length(D);
lyap(zeros(n),D,-C);


Example:

n = 4;
B = rand(n);
C = rand(n);
C*inv(B)-lyap(zeros(n),B,-C)

ans =

1.0e-014 *

-0.2220   -0.2220    0.1221    0.3997
-0.0666   -0.1582   -0.0278   -0.0611
0.0666    0.0555    0.0444    0.5107
0.0222   -0.0222    0.1332    0.7883


This has the benefit of the fact that the code that ships with lyap calls directly into LAPACK routines that should be capable of calling the appropriate sparse solvers/etc.

-
This solution is inelegant, but it is probably the easiest to employ out of the box. I don't have the control systems toolbox installed on a machine that can handle matrices of your size, so I have no idea how practical it might be. No claims are made that this is the best solution. – Emily Oct 29 '12 at 22:01

How about this: rewrite $XD = C$ as $D^T X^T = C^T$. Now solve for $X^T$ column by column.

For example, to get the first column of $X^T$, you would solve $D^T x = c_1$, where $c_1$ is the first column of $C^T$.

You will have to solve many systems of equations (one for each column of $X^T$), but that seems unavoidable, and at least each system is tractable.

If you use a method that first factors $D^T$, then (after the expensive factorization step, which is a one time cost) each system can be solved relatively cheaply using this factorization. (But I'm not sure if a factorization method is appropriate for your large problem.)

-