How to work with quite large matrices?

Hi,

i`m working with the state-space matrices (sparse format) for a 3d FEM model which become quite large (DOF >= 400k).

My System is a first order ODE

with quadratic matrices . In order to perform Model order reduction (Arnoldi algorithm) i need to rewrite my system in a form like

where obviously K and L are obtained by leftdivision of A within the equation. Unfortunately i do encounter big problems with the memory due to the size of the matrices.

I`ve already tried to decompose the matrix A using different decomposition algorithms, but nothing suceed so far. How to deal with such a problem?

Tall arrays can also not be applied, since the matrices are sparse and way to large to be computed as full matrices.

Thanks,

Vierl

[Torsten]

No iterative solver available to solve

A*K = B

for K ?

[Dominik V]

Since B is also a matrix i would need to apply the iterative solvers for any row of the matrix B.

I do not know if the solver configurations are sufficient to handle the problem such that if i decompose B=[b1 b2 ... bn] then the solution is just K=[k1 ... kn]...

If it works that way this could be a possible way of dealing with the problem.

Unfortunately it doesn`t work, since the algorithms do not converge...

[Torsten]

You only need to get the inverse of A, thus solve V*A = I for V. There are iterative methods for this:

https://www.hindawi.com/journals/ijmms/2012/134653/

But I think order reduction can directly be applied to your original system:

http://www.csc.kth.se/utbildning/forskar/avhandlingar/doktor/2005_2006/OlssonHenrik.pdf

Best wishes

Torsten.

[Dominik V]

Thanks a lot. I`ll try this out.

Regards,

Dominik

[Dominik V]

Unfortunetely the matrix has to be inverted, i. e. the linear system has to be solved. The MOR needs the A^(-1)B structure for the moment matching of the transfer function.

The only way is to try some iterative solvers.

Answer 1

[Christine Tobler]

The matrices A\B and A\C are very likely to be dense, even if A, B, C are sparse matrices. Since these are of size 400'000, no matter the method you will run into memory problems, and probably need to use a distributed system and the Parallel Computing Toolbox.

However, are you sure you can't modify the Arnoldi method to directly apply to the generalized system with A, B, and C? This paper, Introduction to Krylov Subspace Methods in Model Order Reduction, seems to say that is possible (see page 4).

[Dominik V]

Well, unfortunately the Krylov-Subspaces for which the transformation matrices are constructed needs always a structure where i need to invert the system either to right or left-hand representation.

At page 5 of your paper you can see that they are using as basic vectors (matrices) the A^(-1)b and A^(-1)*E.

[Torsten]

But in many cases, A is singular. So you say that your model reduction method cannot be applied in this case ?

[Dominik V]

At least the projection-based methods like Arnoldi need this structure. If one uses some kind of open-source software it may be possible to organize the matrices in a suitable way for the reduction. Unfortunately with COMSOL (Livelink for Matlab) i`m only able to extract the matrices with mass term (as mentioned above). So it is not possible to reduce them using e. g. Arnoldi

[Christine Tobler]

Sorry for the late follow-up. Is the problem with computing the decomposition of A that this goes out-of-memory, or that A is singular? Or is it that decomposition A works, but then A\K and A\L can't be stored in-memory?

Some ideas, maybe one of them will help:

• There's no need to compute A\K, storing a decomposition of A and the matrix K will allow the computation of A\(K*x) for each step of building the Krylov subspace.
• Maybe the storage needed for a decomposition of A (LU or LDL, I assume) can be reduced by first applying a reordering to the matrix (see functions amd and dissect).
• If A is singular, is it possible that this is only caused by all-zero rows and columns which could be removed in a pre-processing step?
• If all else fails, you could build a new Krylov subspace vector by replacing the A\v step with an iterative method (this would be a bit painful because there are so many parameters to choose at that point).

[Dominik V]

No problem, thanks for the response.

The problem occurs when i`m trying to decompose A, because the computer runs out of memory. I`ve tried a lot of different decomposition algorithms, but always the same problem.

I also tried to save the matrix as a tall array, which does not work because A is sparse (and full runs out of memory). Also it is now possible to decompose A in blocks and put those together after going to tall arrays, because then the size of the tall arrays decreases because of memory issus ...

Thanks a lot, i will try to follow the suggested solutions. Maybe it will work then.

Best regards,

Dominik