Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

minor fix and a unit test for NNLSSolver #271

Merged
merged 8 commits into from
Mar 5, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 2 additions & 0 deletions .github/workflows/run_tests/action.yml
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@ runs:
mpirun -n 3 --oversubscribe tests/test_StaticSVD
./tests/test_IncrementalSVDBrand
mpirun -n 3 --oversubscribe tests/test_IncrementalSVDBrand
./tests/test_NNLS
mpirun -n 3 --oversubscribe tests/test_NNLS
shell: bash
- name: Basis dataset update test
run: |
Expand Down
1 change: 1 addition & 0 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -275,6 +275,7 @@ if(GTEST_FOUND)
IncrementalSVD
IncrementalSVDBrand
GreedyCustomSampler
NNLS
basis_conversion)
foreach(stem IN LISTS unit_test_stems)
add_executable(test_${stem} unit_tests/test_${stem}.cpp)
Expand Down
30 changes: 23 additions & 7 deletions lib/linalg/NNLS.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -58,15 +58,18 @@ extern "C" {
NNLSSolver::NNLSSolver(double const_tol, int min_nnz, int max_nnz,
int verbosity,
double res_change_termination_tol,
double zero_tol, int n_outer, int n_inner)
double zero_tol, int n_outer, int n_inner,
const NNLS_termination criterion)
dreamer2368 marked this conversation as resolved.
Show resolved Hide resolved
: const_tol_(const_tol), min_nnz_(min_nnz), max_nnz_(max_nnz),
verbosity_(verbosity),
res_change_termination_tol_(res_change_termination_tol),
zero_tol_(zero_tol), n_outer_(n_outer), n_inner_(n_inner),
n_proc_max_for_partial_matrix_(15),
NNLS_qrres_on_(false),
qr_residual_mode_(QRresidualMode::hybrid)
qr_residual_mode_(QRresidualMode::hybrid),
d_criterion(criterion)
{
CAROM_VERIFY((d_criterion == NNLS_termination::L2)
|| (d_criterion == NNLS_termination::LINF));
MPI_Comm_rank(MPI_COMM_WORLD, &d_rank);
MPI_Comm_size(MPI_COMM_WORLD, &d_num_procs);
std::cout << "NNLSSolver init on rank " << d_rank << " out of "
Expand Down Expand Up @@ -147,6 +150,13 @@ void NNLSSolver::solve_parallel_with_scalapack(const Matrix& matTrans,
rhs_halfgap -= rhs_lb;
rhs_halfgap *= 0.5;

double l2norm_threshold;
if (d_criterion == NNLS_termination::L2) {
l2norm_threshold = rhs_halfgap.norm();
if (d_rank == 0 && verbosity_ > 1)
printf("L2 norm threshold: %.5e\n", l2norm_threshold);
}

Vector rhs_avg_glob(rhs_avg);
Vector rhs_halfgap_glob(rhs_halfgap);

Expand All @@ -159,7 +169,7 @@ void NNLSSolver::solve_parallel_with_scalapack(const Matrix& matTrans,
char notrans = 'N';
char layout = 'R';

int n_proc = std::min(n_proc_max_for_partial_matrix_,d_num_procs);
int n_proc = d_num_procs;
int nb = 3; // block column size

// initialize blacs process grid
Expand All @@ -168,8 +178,8 @@ void NNLSSolver::solve_parallel_with_scalapack(const Matrix& matTrans,

int n_dist_loc_max = 0; // maximum number of columns in distributed matrix
if (d_rank < n_proc) {
// m is the maximum number of columns in the global matrix
n_dist_loc_max = ((m/nb + 1)/n_proc + 1)*nb;
// n_tot is the maximum number of columns in the global (not transposed) matrix
n_dist_loc_max = ((n_tot/nb + 1)/n_proc + 1)*nb;
}

std::vector<double> mu_max_array(d_num_procs);
Expand Down Expand Up @@ -254,6 +264,7 @@ void NNLSSolver::solve_parallel_with_scalapack(const Matrix& matTrans,

double mumax_glob;
double rmax;
bool tolerance_met;

for (unsigned int oiter = 0; oiter < n_outer_; ++oiter) {
stalledFlag = 0;
Expand All @@ -265,12 +276,17 @@ void NNLSSolver::solve_parallel_with_scalapack(const Matrix& matTrans,
// This norm of a non-distributed vector is identical on all ranks.
l2_res_hist(oiter) = res_glob.norm();

if (d_criterion == NNLS_termination::LINF)
tolerance_met = (rmax <= const_tol_);
else if (d_criterion == NNLS_termination::L2)
tolerance_met = (l2_res_hist(oiter) <= l2norm_threshold);

if (verbosity_ > 1 && d_rank == 0) {
printf("%d %d %d %d %d %.15e %.15e\n", oiter, n_total_inner_iter,
m, n_tot, n_glob, rmax, l2_res_hist(oiter));
fflush(stdout);
}
if (rmax <= const_tol_ && n_glob >= min_nnz_cap) {
if (tolerance_met && n_glob >= min_nnz_cap) {
if (d_rank == 0 && verbosity_ > 1) {
printf("target tolerance met\n");
fflush(stdout);
Expand Down
21 changes: 14 additions & 7 deletions lib/linalg/NNLS.h
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,16 @@

namespace CAROM {

/**
* @brief Termination criterion for NNLS solver:
* LINF: L-infinity norm (maximum value) of the residual
* L2: L2 norm of the residual
*/
enum class NNLS_termination {
LINF,
L2
};

/**
* \class NNLSSolver
* Class for solving non-negative least-squares problems, cf. T. Chapman et al,
Expand All @@ -34,7 +44,8 @@ class NNLSSolver {
int verbosity=0,
double res_change_termination_tol=1.0e-4,
double zero_tol=1.0e-14, int n_outer=100000,
int n_inner=100000);
int n_inner=100000,
const NNLS_termination criterion=NNLS_termination::LINF);

/**
* Destructor*/
Expand Down Expand Up @@ -115,19 +126,15 @@ class NNLSSolver {
*/
double res_change_termination_tol_;

/**
* @brief Maximum number of processors used in the partial matrix containing
* only the nonzero quadrature points.
*/
int n_proc_max_for_partial_matrix_;

bool normalize_const_;
bool QR_reduce_const_;
bool NNLS_qrres_on_;
QRresidualMode qr_residual_mode_;

int d_num_procs;
int d_rank;

NNLS_termination d_criterion;
};

}
Expand Down
244 changes: 244 additions & 0 deletions unit_tests/test_NNLS.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,244 @@
/******************************************************************************
*
* Copyright (c) 2013-2024, Lawrence Livermore National Security, LLC
* and other libROM project developers. See the top-level COPYRIGHT
* file for details.
*
* SPDX-License-Identifier: (Apache-2.0 OR MIT)
*
*****************************************************************************/

#ifdef CAROM_HAS_GTEST

#include<gtest/gtest.h>
#include "linalg/NNLS.h"
#include <algorithm>
#include <cmath>
#include <cstdio>
#include <cstring> // for memcpy
#include <random>
#include "mpi.h"
#include "utils/mpi_utils.h"

/**
* Simple smoke test to make sure Google Test is properly linked
*/
TEST(GoogleTestFramework, GoogleTestFrameworkFound) {
SUCCEED();
}

TEST(NNLS, solve_with_LINF)
{
int nproc;
int rank;
MPI_Comm_size(MPI_COMM_WORLD, &nproc);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);

const int nrow = 50;
const int ncol = 200;
const int ncol_local = CAROM::split_dimension(ncol);
std::vector<int> row_offset(nproc + 1);
const int total_cols = CAROM::get_global_offsets(ncol_local, row_offset,
MPI_COMM_WORLD);
const double rel_tol = 0.05;
const double nnls_tol = 1.0e-11;
const int true_nnz = 44;

std::default_random_engine generator;
generator.seed(
1234); // fix the seed to keep the same result for different nproc.
std::uniform_real_distribution<> uniform_distribution(0.0, 1.0);
std::normal_distribution<double> normal_distribution(0.0, 1.0);

// distribute from a global matrix to keep the same system for different nproc.
CAROM::Matrix Gt(ncol, nrow, false);
for (int i = 0; i < ncol; i++)
for (int j = 0; j < nrow; j++)
Gt(i, j) = normal_distribution(generator);
Gt.distribute(ncol_local);

CAROM::Vector fom_sol(ncol_local, true);
CAROM::Vector rom_sol(ncol_local, true);
CAROM::Vector rhs(nrow, false);

// distribute from a global matrix to keep the same system for different nproc.
CAROM::Vector fom_sol_serial(ncol, false);
for (int c = 0; c < ncol; c++)
fom_sol_serial(c) = uniform_distribution(generator);
for (int c = 0; c < ncol_local; c++)
fom_sol(c) = fom_sol_serial(row_offset[rank] + c);

Gt.transposeMult(fom_sol, rhs);
rom_sol = 0.0;

CAROM::Vector rhs_lb(rhs);
CAROM::Vector rhs_ub(rhs);

for (int i = 0; i < rhs.dim(); ++i)
{
double delta = rel_tol * abs(rhs(i));
rhs_lb(i) -= delta;
rhs_ub(i) += delta;
}

CAROM::NNLSSolver nnls(nnls_tol, 0, 0, 2);
nnls.solve_parallel_with_scalapack(Gt, rhs_lb, rhs_ub, rom_sol);

int nnz = 0;
for (int i = 0; i < rom_sol.dim(); ++i)
{
if (rom_sol(i) != 0.0)
{
nnz++;
}
}

std::cout << rank << ": Number of nonzeros in NNLS solution: " << nnz
<< ", out of " << rom_sol.dim() << std::endl;

MPI_Allreduce(MPI_IN_PLACE, &nnz, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

if (rank == 0)
std::cout << "Global number of nonzeros in NNLS solution: " << nnz << std::endl;
dreamer2368 marked this conversation as resolved.
Show resolved Hide resolved
EXPECT_EQ(nnz, true_nnz);

// Check residual of NNLS solution
CAROM::Vector res(Gt.numColumns(), false);
Gt.transposeMult(rom_sol, res);

const double normGsol = res.norm();
const double normRHS = rhs.norm();

res -= rhs;
const double relNorm = res.norm() / std::max(normGsol, normRHS);
std::cout << rank << ": relative residual norm for NNLS solution of Gs = Gw: "
<< relNorm << std::endl;

double max_error = 0.0;
for (int k = 0; k < res.dim(); k++)
{
max_error = std::max(max_error, abs(res(k) / rhs(k)));
// printf("rank %d, error(%d): %.3e\n", rank, k, abs(res(k) / rhs(k)));
EXPECT_TRUE(abs(res(k)) < rel_tol * abs(rhs(k)) + nnls_tol);
}
if (rank == 0)
printf("maximum error: %.5e\n", max_error);
}

TEST(NNLS, solve_with_L2)
{
int nproc;
int rank;
MPI_Comm_size(MPI_COMM_WORLD, &nproc);
MPI_Comm_rank(MPI_COMM_WORLD, &rank);

const int nrow = 30;
const int ncol = 100;
const int ncol_local = CAROM::split_dimension(ncol);
std::vector<int> row_offset(nproc + 1);
const int total_cols = CAROM::get_global_offsets(ncol_local, row_offset,
MPI_COMM_WORLD);
const double rel_tol = 0.05;
const double nnls_tol = 1.0e-11;
const int true_nnz = 21;

std::default_random_engine generator;
generator.seed(
1234); // fix the seed to keep the same result for different nproc.
std::uniform_real_distribution<> uniform_distribution(0.0, 1.0);
std::normal_distribution<double> normal_distribution(0.0, 1.0);

// distribute from a global matrix to keep the same system for different nproc.
CAROM::Matrix Gt(ncol, nrow, false);
for (int i = 0; i < ncol; i++)
for (int j = 0; j < nrow; j++)
Gt(i, j) = normal_distribution(generator);
Gt.distribute(ncol_local);

CAROM::Vector fom_sol(ncol_local, true);
CAROM::Vector rom_sol(ncol_local, true);
CAROM::Vector rhs(nrow, false);

// distribute from a global matrix to keep the same system for different nproc.
CAROM::Vector fom_sol_serial(ncol, false);
for (int c = 0; c < ncol; c++)
fom_sol_serial(c) = uniform_distribution(generator);
for (int c = 0; c < ncol_local; c++)
fom_sol(c) = fom_sol_serial(row_offset[rank] + c);

Gt.transposeMult(fom_sol, rhs);
rom_sol = 0.0;

CAROM::Vector rhs_lb(rhs);
CAROM::Vector rhs_ub(rhs);

for (int i = 0; i < rhs.dim(); ++i)
{
double delta = rel_tol * abs(rhs(i));
rhs_lb(i) -= delta;
rhs_ub(i) += delta;
}

CAROM::NNLSSolver nnls(nnls_tol, 0, 0, 2, 1.0e-4, 1.0e-14, 100000, 100000,
CAROM::NNLS_termination::L2);
nnls.solve_parallel_with_scalapack(Gt, rhs_lb, rhs_ub, rom_sol);

int nnz = 0;
for (int i = 0; i < rom_sol.dim(); ++i)
{
if (rom_sol(i) != 0.0)
{
nnz++;
}
}

std::cout << rank << ": Number of nonzeros in NNLS solution: " << nnz
<< ", out of " << rom_sol.dim() << std::endl;

MPI_Allreduce(MPI_IN_PLACE, &nnz, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD);

if (rank == 0)
std::cout << "Global number of nonzeros in NNLS solution: " << nnz << std::endl;
EXPECT_EQ(nnz, true_nnz);

// Check residual of NNLS solution
CAROM::Vector res(Gt.numColumns(), false);
Gt.transposeMult(rom_sol, res);

const double normGsol = res.norm();
const double normRHS = rhs.norm();

res -= rhs;
const double relNorm = res.norm() / std::max(normGsol, normRHS);
std::cout << rank << ": relative residual norm for NNLS solution of Gs = Gw: "
<< relNorm << std::endl;

EXPECT_TRUE(res.norm() < rel_tol * normRHS + nnls_tol);
double max_error = 0.0;
for (int k = 0; k < res.dim(); k++)
{
max_error = std::max(max_error, abs(res(k) / rhs(k)));
}
if (rank == 0)
printf("maximum error: %.5e\n", max_error);
}

int main(int argc, char* argv[])
{
::testing::InitGoogleTest(&argc, argv);
MPI_Init(&argc, &argv);
int result = RUN_ALL_TESTS();
MPI_Finalize();
return result;
}

#else // #ifndef CAROM_HAS_GTEST

int main()
{
std::cout << "libROM was compiled without Google Test support, so unit "
<< "tests have been disabled. To enable unit tests, compile "
<< "libROM with Google Test support." << std::endl;
}

#endif // #endif CAROM_HAS_GTEST
Loading