1

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?

Chris Taylor
  • 28,955
user34790
  • 4,192
  • How big is a "huge" matrix? – Chris Taylor Oct 11 '12 at 11:44
  • 1
    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
  • 2
    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
  • To elaborate slightly on @ChrisTaylor 's comment, the matrix eqn $XD = C$ can be solved for one row of $X$ using one corresponding row of $C$ at a time. There are iterative approximations for $D^{-1}$ but computationally expensive like $O(n^3)$ per iteration, and even if $D^{-1}$ were handed to us, computing $CD^{-1}$ is also $O(mn^2)$ where $m$ is number of rows in $C$ (equiv. in $X$). – hardmath Oct 11 '12 at 13:42
  • @Chris. Yeah I can use that so that I can solve it for one column of X at a time. But If I have to iterate for thousand of times, it's gonna cost me a lot of time, many hours – user34790 Oct 11 '12 at 15:43
  • @user34790 Have you tried block methods? See these notes p. 16 onwards for a description with multiple RHS vectors. – Daryl Oct 11 '12 at 20:26
  • If you have the Control Systems toolbox, you might try the following:

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

     – Emily Oct 29 '12 at 21:51
  • What kind of special structure does $D^T$ have beyond sparsity? Is it symmetric positive definite or anything like that? – littleO Oct 29 '12 at 23:25

3 Answers3

2

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

Rorrik
  • 21
  • 1
  • 2
    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
1

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.

Emily
  • 35,688
  • 6
  • 93
  • 141
  • 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
0

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

littleO
  • 51,938