Skip to content

Handle rank more carefully in sparse GP solves #332

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.

Sign up for GitHub

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

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Handle rank more carefully in sparse GP solves #332

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
4 changes: 3 additions & 1 deletion include/albatross/src/models/sparse_gp.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -315,6 +315,7 @@ class SparseGaussianProcessRegression
}
B.bottomRows(n_new) = A_ldlt.sqrt_solve(K_fu);
const auto B_qr = B.colPivHouseholderQr();
ALBATROSS_ASSERT(B_qr.rank() == B_qr.cols());

// Form:
// y_aug = |R_old P_old^T v_old|
Expand Down Expand Up @@ -379,8 +380,9 @@ class SparseGaussianProcessRegression
compute_internal_components(u, features, targets, &A_ldlt, &K_uu_ldlt,
&K_fu, &y);
const auto B_qr = compute_sigma_qr(K_uu_ldlt, A_ldlt, K_fu);
ALBATROSS_ASSERT(B_qr.rank() == B_qr.cols());

Eigen::VectorXd y_augmented = Eigen::VectorXd::Zero(B_qr.matrixR().rows());
Eigen::VectorXd y_augmented = Eigen::VectorXd::Zero(B_qr.rows());
if (Base::use_async_) {
y_augmented.topRows(y.size()) = A_ldlt.async_sqrt_solve(y);
} else {
Expand Down
8 changes: 6 additions & 2 deletions include/albatross/src/utils/linalg_utils.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,8 +20,10 @@ get_R(const Eigen::ColPivHouseholderQR<Eigen::MatrixXd> &qr) {
// Unfortunately the matrixR() method in Eigen's QR decomposition isn't
// actually the R matrix, it's tall skinny matrix whose lower trapezoid
// contains internal data, only the upper triangular portion is useful
//
// https://eigen.tuxfamily.org/dox/classEigen_1_1ColPivHouseholderQR.html#abe18678c6af88786f2902eeaf2483251
return qr.matrixR()
.topRows(qr.matrixR().cols())
.topLeftCorner(qr.rank(), qr.rank())
.template triangularView<Eigen::Upper>();
}

Expand All @@ -32,7 +34,9 @@ template <typename MatrixType>
inline Eigen::MatrixXd sqrt_solve(const Eigen::MatrixXd &R,
const Eigen::VectorXi &permutation_indices,
const MatrixType &rhs) {

// If `R` is not full-rank, the solve operation below isn't
// well-defined.
ALBATROSS_ASSERT(R.rows() == rhs.rows());
Eigen::MatrixXd sqrt(rhs.rows(), rhs.cols());
for (Eigen::Index i = 0; i < permutation_indices.size(); ++i) {
sqrt.row(i) = rhs.row(permutation_indices.coeff(i));
Expand Down