3
$\begingroup$

recently I have been working on solving some math problems using Fortran. There occurs to me that a linear matrix equation:

$$ AX-XA=B $$

where $A$ and $B$ are known $n\times n$ matrices and $X$ is the one need to be solved. I know this looks like a typical Lyapunov equation. However, in order to solve this equation space friendly, one can use an iterative way to get the numerical result of $X$ when the Kronecker product:

$$ AI-IA^T $$

is non-singular.

Unfortunately, this is not my case. The Kronecker product in my problem is singular. Thus, I cannot use a space friendly iterative way to solve the problem but use a very space consuming method which generates and stores the Kronecker product explicitly!!! Then I used Moore-Penrose pseudoinverse of matrix algorithm to generate the pseudoinverse of this huge matrix and used a library matrix-vector multiplication routine to solve the equation.

For small matrices, this is alright. But when the size of matrices grows (e.g. when n goes up to several hundred) my computer's memory has been completely used up. Any one can help on this so that I can use a space friendly algorithm for such a problem? Thank you very much.

Improve this question
$\endgroup$
2
  • $\begingroup$ As you said, the operator is singular, so there are either infinite or no solutions. Are you fine with any solution, or do you look for one in particular? Which one? $\endgroup$
    – Federico Poloni
    Apr 10, 2014 at 15:47
  • $\begingroup$ @ Federico Poloni To be honest, as long as there is a solution, I can use it directly. So is there a way that generate any solution without using such large memory? $\endgroup$
    – P_E_M_Lee
    Apr 11, 2014 at 3:51

3 Answers 3

Reset to default
0
$\begingroup$

Another idea that might help is mimicking the Bartels-Stewart algorithm for solving general Sylvester equations.

I am assuming that $A$ and $B$ are full $n\times n$ matrices here, so there is no sparsity property to exploit. The resulting algorithm will cost $O(n^3)$.

You need to compute Schur forms $A=QUQ^T$ and $A=RLR^T$, with $U$ upper triangular, $L$ lower triangular, $diag(U)=diag(L)$ and $Q,R$ orthogonal. (You can get the latter from the Schur form of $A^T$.) Then you can transform your problem into the equation $UY-YL=F$, with $Y=Q^TXR$, $F=Q^TBR$. Using the fact that $U$ and $L$ are triangular, you can solve the system one entry of $Y$ at a time with a sort of back-substitution. For instance, comparing the $(n,n)$ entries of the LHS and RHS you get $U_{nn}Y_{nn}-Y_{nn}L_{nn}=F_{nn}$, which has only one unknown and can be solved explicitly. Then you solve for $Y_{n-1,n}$ and so on.

The difference from the "traditional" Bartels-Stewart algorithm is that in your case the equations corresponding to terms on the diagonal will be degenerate, since $U_{nn}=L_{nn}$. So you'll need $F_{ii}=0$ to ensure solvability, and you can set the corresponding elements of $Y$ arbitrarily.

I hope I explained myself -- if you want to learn more on the Bartels-Stewart algorithm, their original article from the 70s is short and readable in my opinion.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Cool. This is worth trying. I will try this out after I finish my current debugging. $\endgroup$
    – P_E_M_Lee
    Apr 13, 2014 at 8:49
  • $\begingroup$ Hi @Federico Poloni, I have tried the lsqr function in Matlab for this particular problem. It turned out that the accuracy is much worse than the method I am currently using (e.g. SVD the huge singular matrix). Moreover, I also plan to try the Schur Decomposition method you provided but I have a new problem: since Uii=Lii, in my case, I am expecting Q(T)BR will give a F that all Fii are zero, which is not the case. Does this mean I need a rotation operation to make all Fii zero? Cheers. $\endgroup$
    – P_E_M_Lee
    Apr 23, 2014 at 8:29
  • $\begingroup$ @P_E_M_Lee I can imagine that the accuracy will be worse in general, but (especially for large dimensions) the method will be faster than a $O(n^6) SVD$; it's a trade-off. There are parameters to require a larger/smaller accuracy goal in lsqr, though, you can try and tweak them. $\endgroup$
    – Federico Poloni
    Apr 23, 2014 at 10:19
  • $\begingroup$ @P_E_M_Lee As for the zeros in the Schur form, the fact that you get non-zeros should mean that the problem is not exactly solvable (i.e., that $vec(B)$ is not in the range of your singular operator). If (1) you know in advance that the problem is solvable and (2) get something significantly larger than the machine precision on the diagonal, there must be something wrong, either in the theory or the implementation. Otherwise I wouldn't worry and I would just truncate $10^{-14}$'s to zeros artificially. $\endgroup$
    – Federico Poloni
    Apr 23, 2014 at 10:22
  • $\begingroup$ Alright then @Federico Poloni, I will try to play with the lsqr for more accuracy. Meanwhile, my boss suggested to regularize the ill-conditioned problem. So that we can transform the singular problem to a normal one and to use an efficient iterative algorithm to calculate the X. Any suggestions on how to regularize the problem? I am still looking for documents related to this topic. $\endgroup$
    – P_E_M_Lee
    Apr 24, 2014 at 6:21
0
$\begingroup$

The solve_sylvester function in scipy will solve this. I don't know whether or not the singularity you mention would affect it.

Improve this answer
$\endgroup$
1
  • $\begingroup$ I have checked a 2x2 example, and solve_sylvester returns an unusable solution with entries of the order of $\mathbf{u}^{-1}$ (machine precision). $\endgroup$
    – Federico Poloni
    Apr 11, 2014 at 12:30
0
$\begingroup$

If a least-squares solution to your singular problem suffices, you might try using the iterative methods in Iterative least-squares solutions of coupled Sylvester equations. These methods should require less memory than your current approach.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Right. I have checked the paper in a short time. It seems the non-singular condition holds for an iterative solution. For singular case, it looks like still has no effective solution? Maybe I am wrong as I did not have time to check the paper carefully. I will read through it more carefully later when I have some time. $\endgroup$
    – P_E_M_Lee
    Apr 11, 2014 at 2:23
  • 1
    $\begingroup$ There are iterative methods for singular systems, e.g., lsqr (stanford.edu/group/SOL/software/lsqr, mathworks.it/it/help/matlab/ref/lsqr.html). $\endgroup$
    – Federico Poloni
    Apr 11, 2014 at 6:20
  • $\begingroup$ @Federico Poloni Tanks for the direction. I will check it carefully later. But I would like to ask, for solving the Ax=b equation of this method, is it necessary to store the large and sparse matrix A explicitly? $\endgroup$
    – P_E_M_Lee
    Apr 11, 2014 at 9:19
  • 1
    $\begingroup$ @P_E_M_Lee No, it's an iterative Krylov-type method. You only need to provide two black-box functions that apply the maps $v\mapsto Mv$ and $w\mapsto M^Tw$, where $M$ is the system matrix. $\endgroup$
    – Federico Poloni
    Apr 11, 2014 at 11:11

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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