search

LLNL / libROM

Model reduction library with an emphasis on large scale parallelism and linear subspace methods
https://www.librom.net
Other
198 stars 36 forks source link

Conventions for the output of SVD #223

Closed dreamer2368 closed 1 year ago

dreamer2368 commented 1 year ago

This is also related to #195 .

I'm currently working on StaticSVD to support the case when the sample size is larger than the system dimension (so the snapshot matrix is transposed). However, there seems to be an inconsistency about the output of SVD, so I wasn't sure which convention to follow.

Is d_basis_right transposed or not?

In StaticSVD::computeSVD and RandomizedSVD::computeSVD, d_basis_right is initialized with the size of transposed matrix (ncolumns, d_num_samples). However, based on the gather_block and gather_transposed_block operations in computeSVD, the output matrix is not supposed to be the transpose, This has never been an issue in the unit tests, since in the tests d_basis_right (in the transposed case, d_basis) is always a square matrix.

Distribution of d_basis/d_basis_right among the processes

In the standard case where the snapshot matrix is not transposed, d_basis is always distributed among the processes, while d_basis_right is always not.

This is not true for RandomizedSVD when the snapshot matrix is transposed. The current implementation of RandomizedSVD::computeSVD simply switches d_basis and d_basis_right, so in this case d_basis_right is distributed among the processes, while d_basis is not.

Based on RandomizedSVD, the convention seems to be that we distribute the 'longer' matrix. In other words, if sample size is smaller than the system dimension, d_basis is distributed, otherwise d_basis_right is distributed. Is this correct? or is this behavior of RandomizedSVD is a bug?

chldkdtn commented 1 year ago

Thank you for looking into this, @dreamer2368 . Do you have a good understanding of how gather_block and gather_transpose_block works? For the first point you made above, it depends if we treat d_basis_right as an orthogonal column matrix or not. If we do, then the size (columns, d_num_samples) is correct. d_basis_right has been used in space-time ROM for transport problems so far, but if I remember correctly, d_basis_right is stored as an orthogonal column matrix, so the size of the matrix is correct, but of course, this is from my memory, which is not reliable. For the second point you made, as I said, d_basis_right mainly used for the space-time ROM where the d_basis_right is served as temporal basis. There, we used to have relatively a small number of time steps, so setting it "distributed" was not making sense. However, if we use it to resolve the issue of he case when the sample size is larger than the system dimension, we definitely need to make it "distributed" as in randomized SVD.

dreamer2368 commented 1 year ago
  1. Concerning the first topic, For a svd operation A = U * S * V^T, I believe at least both StaticSVD and RandomizedSVD return V as d_basis_right, not V^T.
    • The unit test Test_RandomizedSVD checks the result based on this convention.
    • Also RandomizedSVD simply switches d_basis and d_basis_right without any transpose, when the sample size is larger than the system dimension. This means that d_basis_right has to be V, not V^T. In such case the size of d_basis_right should be (d_num_samples, n_columns).
  2. Based on @chldkdtn 's saying, it makes sense to keep d_basis_right not distributed.
    • When the sample size is smaller than the system dimension, it doesn't make sense to distribute d_basis_right.
    • On the other case, we distribute d_basis_right but it will be switched with d_basis, so in the return, it is d_basis that is distributed, not d_basis_right.

Based on this, I think RandomizedSVD::computeSVD will have to be fixed. Does anybody know how it is done for IncrementalSVD, @dylan-copeland and @siuwuncheung ?

chldkdtn commented 1 year ago

For the incremental SVD, it is never assumed that the number of snapshots will be ever greater than the system size because a snapshot enters the incremental SVD in a rank one fashion. Thus, again by my limited memory, I do believe that d_basis_right is never distributed.

By the way, I am confused by what n_columns is and what d_num_samples is. I guess d_num_samples is clear, i.e., the number of snapshot vectors. What is n_columns?

dreamer2368 commented 1 year ago

By the way, I am confused by what n_columns is and what d_num_samples is. I guess d_num_samples is clear, i.e., the number of snapshot vectors. What is n_columns?

n_columns refers to the number of basis modes, I believe. It is determined by either a user-specified number of modes or a fraction in the singular value spectrum. So for A = U * S * V^T operation, V has the size of (d_num_samples, n_columns).

dreamer2368 commented 1 year ago

I checked DMD::constructDMD, and its d_basis_right does follow the convention of (d_num_samples, n_columns). It does not uses any SVD class and directly uses scalapack routines. From a quick look-around, any of DMD classes doesn't seem to use SVD so far.

chldkdtn commented 1 year ago

I briefly checked DMD::constructDMD as well and I think you are right that it does not use SVD class. Instead it calls scalapnack rather directly. For Seung Won's work, incremental SVD in libROM needs to be incorporated to take advantage of the existing implementation in libROM. He will need your help for sure on this.

dreamer2368 commented 1 year ago

This issue is addressed by PR #224 .