14 #include "IncrementalSVDBrand.h"
15 #include "utils/HDFDatabase.h"
24 IncrementalSVDBrand::IncrementalSVDBrand(
26 const std::string& basis_file_name) :
31 d_singular_value_tol(options.singular_value_tol)
33 CAROM_VERIFY(options.singular_value_tol >= 0);
40 if (options.restore_state && d_state_database) {
43 d_state_database->getInteger(
"Up_num_rows", num_rows);
45 d_state_database->getInteger(
"Up_num_cols", num_cols);
46 d_Up =
new Matrix(num_rows, num_cols,
true);
47 d_state_database->getDoubleArray(
"Up",
52 d_state_database->close();
53 delete d_state_database;
60 IncrementalSVDBrand::~IncrementalSVDBrand()
68 if (d_save_state && (!isFirstSample())) {
71 d_state_database->create(d_state_file_name);
74 int num_rows = d_Up->numRows();
75 d_state_database->putInteger(
"Up_num_rows", num_rows);
76 int num_cols = d_Up->numColumns();
77 d_state_database->putInteger(
"Up_num_cols", num_cols);
78 d_state_database->putDoubleArray(
"Up",
90 IncrementalSVDBrand::getSpatialBasis()
94 CAROM_ASSERT(d_basis != 0);
99 IncrementalSVDBrand::getTemporalBasis()
101 updateTemporalBasis();
103 CAROM_ASSERT(d_basis_right != 0);
104 return d_basis_right;
108 IncrementalSVDBrand::buildInitialSVD(
111 CAROM_VERIFY(u != 0);
114 d_S =
new Vector(1,
false);
115 Vector u_vec(u, d_dim,
true);
116 double norm_u = u_vec.norm();
117 d_S->item(0) = norm_u;
120 d_Up =
new Matrix(1, 1,
false);
121 d_Up->item(0, 0) = 1.0;
124 d_U =
new Matrix(d_dim, 1,
true);
125 for (
int i = 0; i < d_dim; ++i) {
126 d_U->item(i, 0) = u[i]/norm_u;
130 if (d_update_right_SV) {
131 d_W =
new Matrix(1, 1,
false);
132 d_W->item(0, 0) = 1.0;
142 IncrementalSVDBrand::buildIncrementalSVD(
143 double* u,
bool add_without_increase)
145 CAROM_VERIFY(u != 0);
149 Vector u_vec(u, d_dim,
true);
150 Vector e_proj(u, d_dim,
true);
151 e_proj -= *(d_U->mult(d_U->transposeMult(e_proj)));
152 e_proj -= *(d_U->mult(d_U->transposeMult(e_proj)));
154 double k = e_proj.inner_product(e_proj);
156 if(d_rank == 0) printf(
"linearly dependent sample!\n");
165 bool linearly_dependent_sample;
166 if ( k < d_linearity_tol ) {
168 std::cout <<
"linearly dependent sample! k = " << k <<
"\n";
169 std::cout <<
"d_linearity_tol = " << d_linearity_tol <<
"\n";
172 linearly_dependent_sample =
true;
173 }
else if ( d_num_samples >= d_max_basis_dimension || add_without_increase ) {
175 linearly_dependent_sample =
true;
183 else if (d_num_samples >= d_total_dim) {
184 linearly_dependent_sample =
true;
187 linearly_dependent_sample =
false;
192 Vector* U_mult_u =
new Vector(d_U->transposeMult(u_vec)->getData(),
195 Vector* l = d_Up->transposeMult(U_mult_u);
203 bool result = svd(Q, A, sigma, W);
214 if ((linearly_dependent_sample && !d_skip_linearly_dependent) ) {
217 if(d_rank == 0) std::cout <<
"adding linearly dependent sample!\n";
218 addLinearlyDependentSample(A, W, sigma);
221 else if (!linearly_dependent_sample) {
225 Vector* j =
new Vector(e_proj.getData(), d_dim,
false);
226 for (
int i = 0; i < d_dim; ++i) {
232 addNewSample(j, A, W, sigma);
247 IncrementalSVDBrand::updateSpatialBasis()
249 d_basis = d_U->mult(d_Up);
252 if ( (d_singular_value_tol != 0.0) &&
253 (d_S->item(d_num_samples-1) < d_singular_value_tol) &&
254 (d_num_samples != 1) ) {
258 "removing a spatial basis corresponding to the small singular value!\n";
261 Matrix* d_basis_new =
new Matrix(d_dim, d_num_samples-1,
262 d_basis->distributed());
263 for (
int row = 0; row < d_dim; ++row) {
264 for (
int col = 0; col < d_num_samples-1; ++col) {
265 d_basis_new->item(row, col) = d_basis->item(row,col);
269 d_basis = d_basis_new;
274 if (fabs(checkOrthogonality(d_basis)) >
275 std::numeric_limits<double>::epsilon()*
static_cast<double>(d_num_samples)) {
276 d_basis->orthogonalize();
282 IncrementalSVDBrand::updateTemporalBasis()
284 delete d_basis_right;
285 d_basis_right =
new Matrix(*d_W);
288 if ( (d_singular_value_tol != 0.0) &&
289 (d_S->item(d_num_samples-1) < d_singular_value_tol) &&
290 (d_num_samples != 1) ) {
294 "removing a temporal basis corresponding to the small singular value!\n";
297 Matrix* d_basis_right_new =
new Matrix(d_num_rows_of_W, d_num_samples-1,
298 d_basis_right->distributed());
299 for (
int row = 0; row < d_num_rows_of_W; ++row) {
300 for (
int col = 0; col < d_num_samples-1; ++col) {
301 d_basis_right_new->item(row, col) = d_basis_right->item(row,col);
304 delete d_basis_right;
305 d_basis_right = d_basis_right_new;
310 if (fabs(checkOrthogonality(d_basis_right)) >
311 std::numeric_limits<double>::epsilon()*d_num_samples) {
312 d_basis_right->orthogonalize();
318 IncrementalSVDBrand::computeBasis()
321 std::cout <<
"d_num_samples = " << d_num_samples <<
"\n";
322 std::cout <<
"d_num_rows_of_W = " << d_num_rows_of_W <<
"\n";
323 std::cout <<
"d_singular_value_tol = " << d_singular_value_tol <<
"\n";
324 std::cout <<
"smallest SV = " << d_S->item(d_num_samples-1) <<
"\n";
325 if (d_num_samples > 1) {
326 std::cout <<
"next smallest SV = " << d_S->item(d_num_samples-2) <<
"\n";
330 updateSpatialBasis();
331 if (d_update_right_SV)
333 updateTemporalBasis();
337 if ( (d_singular_value_tol != 0.0) &&
338 (d_S->item(d_num_samples-1) < d_singular_value_tol) &&
339 (d_num_samples != 1) ) {
346 IncrementalSVDBrand::addLinearlyDependentSample(
351 CAROM_VERIFY(A != 0);
352 CAROM_VERIFY(sigma != 0);
356 Matrix Amod(d_num_samples, d_num_samples,
false);
357 for (
int row = 0; row < d_num_samples; ++row) {
358 for (
int col = 0; col < d_num_samples; ++col) {
359 Amod.item(row, col) = A->item(row, col);
362 d_S->item(col) = sigma->item(row, col);
368 Matrix* Up_times_Amod = d_Up->mult(Amod);
370 d_Up = Up_times_Amod;
373 if (d_update_right_SV) {
379 new_d_W =
new Matrix(d_num_rows_of_W+1, d_num_samples,
false);
380 for (
int row = 0; row < d_num_rows_of_W; ++row) {
381 for (
int col = 0; col < d_num_samples; ++col) {
382 double new_d_W_entry = 0.0;
383 for (
int entry = 0; entry < d_num_samples; ++entry) {
384 new_d_W_entry += d_W->item(row, entry)*W->item(entry, col);
386 new_d_W->item(row, col) = new_d_W_entry;
389 for (
int col = 0; col < d_num_samples; ++col) {
390 new_d_W->item(d_num_rows_of_W, col) = W->item(d_num_samples, col);
400 IncrementalSVDBrand::addNewSample(
406 CAROM_VERIFY(j != 0);
407 CAROM_VERIFY(A != 0);
408 CAROM_VERIFY(sigma != 0);
411 Matrix* newU =
new Matrix(d_dim, d_num_samples+1,
true);
412 for (
int row = 0; row < d_dim; ++row) {
413 for (
int col = 0; col < d_num_samples; ++col) {
414 newU->item(row, col) = d_U->item(row, col);
416 newU->item(row, d_num_samples) = j->item(row);
422 if (d_update_right_SV) {
428 new_d_W =
new Matrix(d_num_rows_of_W+1, d_num_samples+1,
false);
429 for (
int row = 0; row < d_num_rows_of_W; ++row) {
430 for (
int col = 0; col < d_num_samples+1; ++col) {
431 double new_d_W_entry = 0.0;
432 for (
int entry = 0; entry < d_num_samples; ++entry) {
433 new_d_W_entry += d_W->item(row, entry)*W->item(entry, col);
435 new_d_W->item(row, col) = new_d_W_entry;
438 for (
int col = 0; col < d_num_samples+1; ++col) {
439 new_d_W->item(d_num_rows_of_W, col) = W->item(d_num_samples, col);
450 Matrix* new_d_Up =
new Matrix(d_num_samples+1, d_num_samples+1,
false);
451 for (
int row = 0; row < d_num_samples; ++row) {
452 for (
int col = 0; col < d_num_samples+1; ++col) {
453 double new_d_Up_entry = 0.0;
454 for (
int entry = 0; entry < d_num_samples; ++entry) {
455 new_d_Up_entry += d_Up->item(row, entry)*A->item(entry, col);
457 new_d_Up->item(row, col) = new_d_Up_entry;
460 for (
int col = 0; col < d_num_samples+1; ++col) {
461 new_d_Up->item(d_num_samples, col) = A->item(d_num_samples, col);
468 int num_dim = std::min(sigma->numRows(), sigma->numColumns());
469 d_S =
new Vector(num_dim,
false);
470 for (
int i = 0; i < num_dim; i++) {
471 d_S->item(i) = sigma->item(i,i);
480 if (d_num_samples > d_total_dim) {
481 max_U_dim = d_num_samples;
484 max_U_dim = d_total_dim;
486 if (fabs(checkOrthogonality(d_Up)) >
487 std::numeric_limits<double>::epsilon()*
static_cast<double>(max_U_dim)) {
488 d_Up->orthogonalize();
490 if (fabs(checkOrthogonality(d_U)) >
491 std::numeric_limits<double>::epsilon()*
static_cast<double>(max_U_dim)) {
492 d_U->orthogonalize();
495 if(d_update_right_SV)
497 if (fabs(checkOrthogonality(d_W)) >
498 std::numeric_limits<double>::epsilon()*d_num_samples) {
499 d_W->orthogonalize();