LinearOperator for fully dense and non-square matrix

[Tao Jin]

Dear all,

I am developing a quasi Newton solver in deal.ii. During each iteration, I need to solve a linear system B * x = p, where the matrix B has the following form:

B = B_0 + W * M * W^T.

B_0 is an n by n sparse matrix (n is large > 100k);

M is a m by m full matrix (m is small m = 20 for example);

W is a n by m full matrix;

Since W * M * W^T will be an n by n fully dense matrix, I cannot afford to explicitly store it. 

Does the deal.ii community have any suggestion what is the most efficient way to solve this linear system? I am thinking to use LinearOperator to define and operate on W * M * W^T, but I am not sure whether it supports full matrix and how efficient it would be.

Best,

Tao

[Wolfgang Bangerth]

On 4/18/24 12:06, Tao Jin wrote:
>
> I am developing a quasi Newton solver in deal.ii. During each iteration, I
> need to solve a linear system B * x = p, where the matrix B has the following
> form:
> B = B_0 + W * M * W^T.
> B_0 is an n by n sparse matrix (n is large > 100k);

> M is a m by m *full matrix* (m is small m = 20 for example);

> W is a n by m full matrix;
> Since W * M * W^T will be an n by n fully dense matrix, I cannot afford to
> explicitly store it.
>
> Does the deal.ii community have any suggestion what is the most efficient way
> to solve this linear system? I am thinking to use LinearOperator to define and
> operate on W * M * W^T, but I am not sure whether it supports full matrix and
> how efficient it would be.

Tao:
I'm not an expert in the LinearOperator framework, but operators such as yours
appear, for example, in step-20:
https://dealii.org/developer/doxygen/deal.II/step_20.html#step_20-Implementationoflinearsolversandpreconditioners
There, opS has exactly the form you describe, and it does not actually compute
the entries of the S matrix, but only represents the *action* as three
matrix-vector products.

That only leaves the question of whether you can use linear_operator(M) on
full matrices M. I don't know the answer to that, but what happens if you try?

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/

[Tao Jin]

Hi Wolfgang,

Thank you so much for your reply. I think my problem has two potential issues:

1. W * M * W^T will be an n by n fully dense matrix with n > 100k. Even if it is valid to use LinearOperator to perform W * M * W^T, will memory required to store this operator be an issue?

2. I still have to solve B x= (B_0 + W * M * W^T) x = p. Whether linear solvers such as CG will be efficient enough is also an issue since W * M * W^T is fully dense.

Anyhow, I will do some experiments and let the community know.

Best,

Tao

[Wolfgang Bangerth]

On 4/19/24 19:42, Tao Jin wrote:
>
> Thank you so much for your reply. I think my problem has two potential issues:
> 1. W * M * W^T will be an n by n fully dense matrix with n > 100k. Even if it
> is valid to use LinearOperator to perform W * M * W^T, will memory required to
> store this operator be an issue?
> 2. I still have to solve B x= (B_0 + W * M * W^T) x = p. Whether linear
> solvers such as CG will be efficient enough is also an issue since W * M * W^T
> is fully dense.

Tao,
like I mentioned, a line like
S = linear_operator(W) * linear_operator(M) * transpose_operator(W)
does not actually compute the product of these matrices. The documentation of
step-20 and step-22 goes to great length in explaining this issue. All the
object S provides is the ability to multiple S by a vector, for which you only
need matrix-vector products because
S x
= (W M W^T) x
= W (M (W^T x))

The product of matrices is never computed because it is never needed.

[Tao Jin]

Dear Wolfgang,

What you said makes sense. 

Now that my LinearOperator is defined as B_0 + W * M * W^T, I only need to find an appropriate iterative solver that can solve the inverse of this operator (the underlying matrix is fully dense) efficiently.

Best,

Tao

[Wolfgang Bangerth]

On 4/19/24 21:54, Tao Jin wrote:
>
> Now that my LinearOperator is defined as B_0 + W * M * W^T, I only need to
> find an appropriate iterative solver that can solve the inverse of this
> operator (the underlying matrix is fully dense) efficiently.

It is entirely unimportant that the matrix is dense. For iterative solvers,
this is of no consequence at all. What matters is whether, for example, it is
symmetric and/or positive definite. The other things that matters is whether
you can construct a good preconditioner for it.

You really should read through step-20 and step-22 to see how all of this
works in practice.

[Tao Jin]

Dear Wolfgang,

Thanks for the advice! 

Best,

Tao

[Tao Jin]

Dear Wolfgang,

I am glad to report that the LinearOperator modules also work for fully dense and non-rectangular matrices (FullMatrix<double>). 

I can now use the inverse operator to solve the linear system:

B = B_0 + W * M * W^T.

B_0 is an n by n sparse symmetric positive-definite matrix (n is large > 100k);

M is a m by m full matrix (m is small m = 20 for example);

W is a n by m full matrix;

Thanks again for your advice.

Best,

Tao