chldkdtn
commented
1 year ago Thank you @dreamer2368 for identifying this failing case. @siuwuncheung can you verify if this is the reason that the tests for your PR fail as well?
Closed dreamer2368 closed 1 year ago
chldkdtn
commented
1 year ago Thank you @dreamer2368 for identifying this failing case. @siuwuncheung can you verify if this is the reason that the tests for your PR fail as well?
siuwuncheung
commented
1 year ago Thank you @dreamer2368 for identifying this failing case. @siuwuncheung can you verify if this is the reason that the tests for your PR fail as well?
The failed test in my PR is a unit test with DMD, so it should be two different issues. I will fix the one in mine today.
[ FAILED ] DMDTest.Test_DMD (43 ms) [----------] 1 test from DMDTest (43 ms total)
[----------] Global test environment tear-down [==========] 2 tests from 2 test suites ran. (43 ms total) [ PASSED ] 1 test. [ FAILED ] 1 test, listed below: [ FAILED ] DMDTest.Test_DMD
siuwuncheung
commented
1 year ago @dreamer2368 I think this does not happen in the CI workflow in GH PRs. Can you please verify?
dreamer2368
commented
1 year ago No it does not happen in the CI workflow. But it can (and does) happen in other machines. That is the main issue, since the test result is susceptible to round-off errors, causing false-negative.
I actually checked that qr decomposition in RandomizedSVD::computeSVD is indeed susceptible to round-off errors. (this only happens with d_debug_algorithm, so only for test.)
Following minimal example is taken from the first test of Test_RandomizedSVD. We see that O(1e-15) perturbation can completely change 2nd/3rd column of Q, due to the rank deficiency caused by rand_mat. To fix this issue, I suggest to change the 'testing' rand_mat and update the answer accordingly.
#include <fstream>
#include <iostream>
#include "librom.h"
// uniform random number in [0., 1.]
double UniformRandom();
int main(int argc, char *argv[])
{
MPI_Init(&argc, &argv);
printf("Example from RandomizedSVDTest, Test_RandomizedSVD.\n\n");
double* sample1 = new double[5] {0.5377, 1.8339, -2.2588, 0.8622, 0.3188};
double* sample2 = new double[5] {-1.3077, -0.4336, 0.3426, 3.5784, 2.7694};
double* sample3 = new double[5] {-1.3499, 3.0349, 0.7254, -0.0631, 0.7147};
CAROM::Matrix *snapshot_matrix = new CAROM::Matrix(5, 3, true);
for (int i = 0; i < snapshot_matrix->numRows(); i++) {
snapshot_matrix->item(i, 0) = sample1[i];
snapshot_matrix->item(i, 1) = sample2[i];
snapshot_matrix->item(i, 2) = sample3[i];
}
printf("Snapshot matrix from example.\n");
for (int i = 0; i < snapshot_matrix->numRows(); i++) {
for (int j = 0; j < snapshot_matrix->numColumns(); j++) {
printf("%.15E\t", snapshot_matrix->item(i, j));
}
printf("\n");
}
printf("\n");
CAROM::Matrix *rand_mat = new CAROM::Matrix(3, 3, false);
for (int i = 0; i < snapshot_matrix->numColumns(); i++) {
for (int j = 0; j < 3; j++) {
rand_mat->item(i, j) = 1;
}
}
printf("rand_mat (%d x %d) in RandomizedSVD::computeSVD (d_debug_algorithm).\n", rand_mat->numRows(), rand_mat->numColumns());
for (int i = 0; i < rand_mat->numRows(); i++) {
for (int j = 0; j < rand_mat->numColumns(); j++) {
printf("%.15E\t", rand_mat->item(i, j));
}
printf("\n");
}
printf("\n");
CAROM::Matrix* rand_proj = snapshot_matrix->mult(rand_mat);
printf("rand_proj (%d x %d) in RandomizedSVD::computeSVD.\n", rand_proj->numRows(), rand_proj->numColumns());
for (int i = 0; i < rand_proj->numRows(); i++) {
for (int j = 0; j < rand_proj->numColumns(); j++) {
printf("%.15E\t", rand_proj->item(i, j));
}
printf("\n");
}
printf("\n");
CAROM::Matrix* Q = rand_proj->qr_factorize();
printf("Q (%d x %d) in RandomizedSVD::computeSVD.\n", Q->numRows(), Q->numColumns());
for (int i = 0; i < Q->numRows(); i++) {
for (int j = 0; j < Q->numColumns(); j++) {
printf("%.15E\t", Q->item(i, j));
}
printf("\n");
}
printf("\n");
// Now perturb the snapshot matrix within machine precision.
for (int i = 0; i < snapshot_matrix->numRows(); i++) {
for (int j = 0; j < snapshot_matrix->numColumns(); j++) {
snapshot_matrix->item(i, j) *= 1.0 + 1.0e-15 * (-1.0 + 2.0 * UniformRandom());
}
}
printf("Perturbed snapshot matrix.\n");
for (int i = 0; i < snapshot_matrix->numRows(); i++) {
for (int j = 0; j < snapshot_matrix->numColumns(); j++) {
printf("%.15E\t", snapshot_matrix->item(i, j));
}
printf("\n");
}
printf("\n");
rand_proj = snapshot_matrix->mult(rand_mat);
Q = rand_proj->qr_factorize();
printf("After perturbation, 2nd/3rd columns of Q totally changed.\n");
for (int i = 0; i < Q->numRows(); i++) {
for (int j = 0; j < Q->numColumns(); j++) {
printf("%.15E\t", Q->item(i, j));
}
printf("\n");
}
printf("\n");
printf("Both Q are legitimate QR decomposition of rand_proj within 'machine precision'.\n");
MPI_Finalize();
return 0;
}
double UniformRandom()
{
return static_cast<double> (rand()) / static_cast<double> (RAND_MAX);
}The resulting output is:
Example from RandomizedSVDTest, Test_RandomizedSVD.
Snapshot matrix from example.
5.377000000000000E-01 -1.307700000000000E+00 -1.349900000000000E+00
1.833900000000000E+00 -4.336000000000000E-01 3.034900000000000E+00
-2.258800000000000E+00 3.426000000000000E-01 7.254000000000000E-01
8.622000000000000E-01 3.578400000000000E+00 -6.310000000000000E-02
3.188000000000000E-01 2.769400000000000E+00 7.147000000000000E-01
rand_mat (3 x 3) in RandomizedSVD::computeSVD (d_debug_algorithm).
1.000000000000000E+00 1.000000000000000E+00 1.000000000000000E+00
1.000000000000000E+00 1.000000000000000E+00 1.000000000000000E+00
1.000000000000000E+00 1.000000000000000E+00 1.000000000000000E+00
rand_proj (5 x 3) in RandomizedSVD::computeSVD.
-2.119900000000000E+00 -2.119900000000000E+00 -2.119900000000000E+00
4.435200000000000E+00 4.435200000000000E+00 4.435200000000000E+00
-1.190800000000000E+00 -1.190800000000000E+00 -1.190800000000000E+00
4.377500000000000E+00 4.377500000000000E+00 4.377500000000000E+00
3.802900000000000E+00 3.802900000000000E+00 3.802900000000000E+00
Q (5 x 3) in RandomizedSVD::computeSVD.
-2.755033519022854E-01 -9.077467635817615E-01 3.135274050896839E-01
5.764009936114987E-01 -4.177980391731229E-01 -7.001788974619734E-01
-1.547570128049631E-01 3.758194159204378E-02 -1.295324450643630E-01
5.689022703675899E-01 -5.221036081088777E-03 5.208590570351792E-01
4.942269432280772E-01 -9.747401413574118E-04 3.512510016878761E-01
Perturbed snapshot matrix.
5.377000000000003E-01 -1.307700000000000E+00 -1.349900000000001E+00
1.833900000000001E+00 -4.336000000000004E-01 3.034899999999998E+00
-2.258799999999999E+00 3.426000000000002E-01 7.253999999999997E-01
8.622000000000000E-01 3.578400000000000E+00 -6.310000000000002E-02
3.187999999999999E-01 2.769400000000000E+00 7.147000000000007E-01
After perturbation, 2nd/3rd columns of Q totally changed.
-2.755033519022854E-01 -9.517690886536551E-01 6.753268486376426E-02
5.764009936114984E-01 -7.711547662313500E-02 7.454636540798325E-01
-1.547570128049630E-01 1.434853849974411E-01 -1.711447593679669E-01
5.689022703675899E-01 -2.064584708791187E-01 -2.607804146307769E-01
4.942269432280773E-01 -1.580368603932945E-01 -5.851720461598588E-01
Both Q are legitimate QR decomposition of rand_proj within 'machine precision'.@dreamer2368 If you can submit a PR with this fix, that would be wonderful.
siuwuncheung
commented
1 year ago Sorry for getting back this late. Thank you so much for setting up the example to make things clear, @dreamer2368. I agree that we need to change the rand_mat. Normally, it is a (Gaussian) random matrix, consisting of linear independent random columns. In the debug mode, we set it to be all ones, which break the linear independence and cause numerical instability issues. We need to be careful about how many tests, besides this one, will be affected by changing the rand_mat, before we know what to change.
If I understand correctly, we always assume the matrix being randomly projected has more or equal rows than columns. In this case, it will be safe to replace rand_mat->item(i, j) = 1; by rand_mat->item(i, j) = (i == j);
dreamer2368
commented
1 year ago Resolved by PR #211 .
The unit test
test_RandomizedSVDfails. This is observed from both docker container librom_env and quartz.This needs more investigations, but the issue seems to be a peculiar setup within the test. For debug/test cases,
rand_matinRandomizedSVD::computeSVDis a matrix of 1. As a result, the projected snapshot matrixrand_projbecomes 1-rank matrix. The resultingQmatrix from qr factorization can have arbitrary orthonormal column basis except one column. The result of SVD will depend on these arbitrary column basis ofQ, although the answer is fixed to one case.If this is the case, then the test
rand_matand the corresponding answer should be adjusted so that it may have some randomness (probably under the fixed random seed).Following output is the result from
test_RandomizedSVD: