Skip to content

Commit bc30c41

Browse files
AlePaluAlePalu
authored
Merge pull request #15 from marcogalliani/stable
Updates c48e9c6, introduced isolsted RSI, RBKI, NysRSI, NysRBKI for random SVD approximation, added related random utils (BCGS, GaussianMayrix)
2 parents af8d503 + 61930c8 commit bc30c41

File tree

6 files changed

+547
-101
lines changed

6 files changed

+547
-101
lines changed

‎fdaPDE/linear_algebra.h

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -29,4 +29,8 @@
2929
#include "linear_algebra/fspai.h"
3030
#include "linear_algebra/lumping.h"
3131

32+
#include "fdaPDE/linear_algebra/randomized_algorithms/rsvd.h"
33+
#include "fdaPDE/linear_algebra/randomized_algorithms/revd.h"
34+
#include "fdaPDE/linear_algebra/randomized_algorithms/nystrom.h"
35+
3236
#endif // __FDAPDE_LINEAR_ALGEBRA_MODULE_H__

‎fdaPDE/linear_algebra/randomized_algorithms/nystrom.h

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
//
2+
// Created by Marco Galliani on 03/11/24.
3+
//
4+
5+
#ifndef NYSTROM_APPROX_H
6+
#define NYSTROM_APPROX_H
7+
8+
#include <memory>
9+
#include <unordered_set>
10+
11+
namespace fdapde{
12+
namespace core{
13+
14+
template<typename MatrixType>
15+
class NysApproxStrategy{
16+
protected:
17+
unsigned int seed_=fdapde::random_seed;
18+
double tol_=1e-3;
19+
//storage of the decomposition
20+
DMatrix<double> F_;
21+
double shift_;
22+
public:
23+
NysApproxStrategy()=default;
24+
NysApproxStrategy(unsigned int seed , double tol) : seed_(seed), tol_(tol){}
25+
virtual void compute(const MatrixType &A, int block_sz, int max_iter) = 0;
26+
//setters
27+
void setTol(double tol){ tol_=tol;}
28+
void setSeed(unsigned int seed){ seed_=seed;}
29+
//getters
30+
int rank() const{ return F_.cols();}
31+
DMatrix<double> factor() const{ return F_;}
32+
double shift() const{ return shift_;}
33+
//destructor
34+
virtual ~NysApproxStrategy() = default;
35+
};
36+
37+
template<typename MatrixType>
38+
class RPChol : public NysApproxStrategy<MatrixType>{
39+
public:
40+
RPChol()=default;
41+
RPChol(unsigned int seed, double tol) : NysApproxStrategy<MatrixType>(seed,tol){}
42+
void compute(const MatrixType &A, int block_sz, int max_iter) override{
43+
//params init
44+
max_iter = std::min(max_iter,(int)std::ceil((double)A.cols()/(double)block_sz));
45+
std::mt19937 rng{this->seed_};
46+
this->shift_ = std::numeric_limits<double>::epsilon()*A.trace();
47+
//factor init
48+
std::vector<int> ind_cols_A(A.cols());
49+
std::iota(ind_cols_A.begin(),ind_cols_A.end(),0);
50+
DVector<double> diag_res = A.diagonal();
51+
this->F_ = DMatrix<double>::Zero(A.rows(),max_iter*block_sz);
52+
//error
53+
double norm_A=A.norm(), reconstruction_err = norm_A;
54+
//iterations
55+
int n_cols_F = 0;
56+
while(n_cols_F<max_iter*block_sz && reconstruction_err>this->tol_*norm_A){
57+
//sampling
58+
std::discrete_distribution<int> sampling_distr(diag_res.begin(),diag_res.end());
59+
std::unordered_set<int> sampled_pivots(block_sz);
60+
for(int j=0; (int)sampled_pivots.size()<block_sz && j<2*block_sz; j++){
61+
sampled_pivots.insert(sampling_distr(rng));
62+
}
63+
std::vector<int> pivot_set(sampled_pivots.begin(),sampled_pivots.end()); //converting for Eigen slicing
64+
//building the F factor
65+
DMatrix<double> G = A(Eigen::all,pivot_set);
66+
G = G - this->F_.leftCols(n_cols_F) * this->F_(pivot_set,Eigen::all).leftCols(n_cols_F).transpose();
67+
double shift = std::numeric_limits<double>::epsilon()*G(pivot_set,Eigen::all).trace();
68+
Eigen::LLT<DMatrix<double>> chol(G(pivot_set,Eigen::all) + shift*DMatrix<double>::Identity(pivot_set.size(),pivot_set.size()));
69+
DMatrix<double> T = chol.matrixU().solve<Eigen::OnTheRight>(G);
70+
this->F_.middleCols(n_cols_F,pivot_set.size()) = T;
71+
//update the sampling distribution
72+
diag_res = (diag_res - T.rowwise().squaredNorm()).array().max(0);
73+
//update the error
74+
n_cols_F += pivot_set.size();
75+
reconstruction_err = (A-this->F_.leftCols(n_cols_F)*this->F_.leftCols(n_cols_F).transpose()).norm();
76+
}
77+
this->F_ = this->F_.leftCols(n_cols_F);
78+
return;
79+
}
80+
};
81+
82+
template<typename MatrixType>
83+
class NystromApproximation{
84+
private:
85+
std::unique_ptr<NysApproxStrategy<MatrixType>> nys_strategy_;
86+
public:
87+
explicit NystromApproximation(std::unique_ptr<NysApproxStrategy<MatrixType>> &&strategy=std::make_unique<RPChol<DMatrix<double>>>()): nys_strategy_(std::move(strategy)){}
88+
void compute(const MatrixType &A, int block_sz, int max_iter=1e3){
89+
nys_strategy_->compute(A,block_sz, max_iter);
90+
return;
91+
}
92+
//setters
93+
void setTol(double tol){ nys_strategy_->setTol(tol);}
94+
void setSeed(unsigned int seed){ nys_strategy_->setSeed(seed);}
95+
//getters
96+
int rank() const{ return nys_strategy_->rank();}
97+
DMatrix<double> factor() const{ return nys_strategy_->factor();}
98+
double shift() const{ return nys_strategy_->shift();}
99+
};
100+
101+
}//core
102+
}//fdpade
103+
104+
#endif //NYSTROM_APPROX_H

‎fdaPDE/linear_algebra/randomized_algorithms/revd.h

Lines changed: 162 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,162 @@
1+
//
2+
// Created by Marco Galliani on 30/10/24.
3+
//
4+
5+
#ifndef REVD_H
6+
#define REVD_H
7+
8+
#include <utility>
9+
#include <memory>
10+
#include <tuple>
11+
#include <limits>
12+
#include <type_traits>
13+
14+
namespace fdapde{
15+
namespace core{
16+
//Interface for the Approximation strategy
17+
template<typename MatrixType>
18+
class REVDStrategy{
19+
protected:
20+
unsigned int seed_=fdapde::random_seed;
21+
double tol_=1e-3;
22+
//storage of the decomposition
23+
DMatrix<double> U_;
24+
DVector<double> Lambda_;
25+
public:
26+
REVDStrategy()=default;
27+
REVDStrategy(unsigned int seed, double tol) : seed_(seed), tol_(tol){}
28+
virtual void compute(const MatrixType &A, int rank, int max_iter) = 0;
29+
//setters
30+
void setTol(double tol){ tol_=tol;}
31+
void setSeed(unsigned int seed){ seed_=seed;}
32+
//getters
33+
int rank() const{ return Lambda_.size();}
34+
DMatrix<double> matrixU() const{ return U_;}
35+
DVector<double> eigenValues() const{ return Lambda_;}
36+
//destructor
37+
virtual ~REVDStrategy() = default;
38+
};
39+
40+
template<typename MatrixType>
41+
class NysRSI : public REVDStrategy<MatrixType>{
42+
public:
43+
NysRSI()=default;
44+
NysRSI(unsigned int seed, double tol) : REVDStrategy<MatrixType>(seed,tol){}
45+
void compute(const MatrixType &A, int rank, int max_iter) override{
46+
//params init
47+
int max_rank = A.rows(); //equal to A.cols()
48+
int block_sz = std::min(2*rank,max_rank); //default setting
49+
max_iter = std::min(max_iter, max_rank);
50+
double shift = A.trace()*std::numeric_limits<double>::epsilon();
51+
//factor init
52+
DMatrix<double> Y = fdapde::internals::GaussianMatrix(A.rows(), block_sz, this->seed_);
53+
DMatrix<double> X;
54+
DMatrix<double> F;
55+
Eigen::HouseholderQR<DMatrix<double>> qr;
56+
//error
57+
Eigen::JacobiSVD<DMatrix<double>> svd;
58+
DMatrix<double> E;
59+
double norm_A = A.norm(), res_err = norm_A;
60+
//iterations
61+
for(int i=0; res_err > this->tol_*norm_A && i<max_iter; ++i) {
62+
qr.compute(Y);
63+
X = qr.householderQ() * DMatrix<double>::Identity(A.rows(),block_sz);
64+
Y = A*X;
65+
//construct the factor
66+
Y += shift*DMatrix<double>::Identity(Y.rows(),Y.cols());
67+
Eigen::LLT<DMatrix<double>> chol(X.transpose()*Y);
68+
F = chol.matrixU().solve<Eigen::OnTheRight>(Y);
69+
//update the error
70+
svd.compute(F,Eigen::ComputeThinU | Eigen::ComputeThinV);
71+
E = A*svd.matrixU().leftCols(rank) - svd.matrixU().leftCols(rank)*(svd.singularValues().head(rank).array().pow(2)-shift).matrix().asDiagonal();
72+
res_err = E.colwise().template lpNorm<2>().maxCoeff();
73+
}
74+
this->U_ = svd.matrixU().leftCols(rank);
75+
this->Lambda_ = (svd.singularValues().head(rank).array().pow(2)-shift).matrix();
76+
return;
77+
}
78+
};
79+
80+
template<typename MatrixType>
81+
class NysRBKI : public REVDStrategy<MatrixType>{
82+
public:
83+
NysRBKI()=default;
84+
NysRBKI(unsigned int seed, double tol) : REVDStrategy<MatrixType>(seed,tol){}
85+
void compute(const MatrixType &A, int rank, int max_iter) override{
86+
//params init
87+
int max_rank = A.rows(); //equal to A.cols()
88+
int block_sz; //default setting
89+
if(A.rows()<1000){
90+
block_sz = 1;
91+
}else{
92+
block_sz = 10;
93+
}
94+
max_iter = std::min(max_iter,max_rank/block_sz-1);
95+
double shift = A.trace()*std::numeric_limits<double>::epsilon();
96+
//factor init
97+
DMatrix<double> X,Y,S,F;
98+
X.resize(A.rows(),max_rank); Y.resize(A.rows(),max_rank);
99+
S = DMatrix<double>::Zero(max_rank,max_rank);
100+
Eigen::HouseholderQR<DMatrix<double>> qr(fdapde::internals::GaussianMatrix(A.rows(),block_sz,this->seed_));
101+
X.leftCols(block_sz) = qr.householderQ()*DMatrix<double>::Identity(A.rows(),block_sz);
102+
Y.leftCols(block_sz) = A*X.leftCols(block_sz);
103+
//error
104+
Eigen::JacobiSVD<DMatrix<double>> svd;
105+
DMatrix<double> E;
106+
double norm_A=A.norm(), res_err=norm_A;
107+
//iterations
108+
int n_cols_X = block_sz;
109+
for(int i=0; i<max_iter && res_err>this->tol_*norm_A;i++,n_cols_X+=block_sz){
110+
X.middleCols((i+1)*block_sz,block_sz) = Y.middleCols(i*block_sz,block_sz) + shift*X.middleCols(i*block_sz,block_sz);
111+
//blocked column
112+
DMatrix<double> new_col = DMatrix<double>::Zero(X.rows(),(i+1)*block_sz);
113+
new_col.middleCols(std::max(i-1,0)*block_sz,block_sz) = X.middleCols(std::max(i-1,0)*block_sz,block_sz);
114+
new_col.middleCols(i*block_sz,block_sz) = X.middleCols(i*block_sz,block_sz);
115+
new_col = new_col.transpose()*X.middleCols((i+1)*block_sz,block_sz);
116+
//orthogonalisation
117+
auto new_block_qr = fdapde::internals::BCGS_plus(X.leftCols((i+1)*block_sz),X.middleCols((i+1)*block_sz,block_sz));
118+
X.middleCols((i+1)*block_sz,block_sz) = new_block_qr.first;
119+
//cholesky
120+
S.block(0,i*block_sz,(i+1)*block_sz,block_sz) = new_col;
121+
Eigen::LLT<DMatrix<double>> chol(S.block(0,0,(i+1)*block_sz,(i+1)*block_sz));
122+
S.block((i+1)*block_sz,i*block_sz,block_sz,block_sz) = new_block_qr.second;
123+
F = chol.matrixU().solve<Eigen::OnTheRight>(S.block(0,0,(i+2)*block_sz,(i+1)*block_sz));
124+
//update Y
125+
Y.middleCols((i+1)*block_sz,block_sz) = A*X.middleCols((i+1)*block_sz,block_sz);
126+
//update the error
127+
svd.compute(F, Eigen::ComputeThinU | Eigen::ComputeThinV);
128+
E = Y.leftCols((i+2)*block_sz)*svd.matrixU().leftCols(std::min(rank,(i+1)*block_sz)) - X.leftCols((i+2)*block_sz)*svd.matrixU().leftCols(std::min(rank,(i+1)*block_sz))*(svd.singularValues().head(std::min(rank,(i+1)*block_sz)).array().pow(2)-shift).matrix().asDiagonal();
129+
res_err = E.colwise().template lpNorm<2>().maxCoeff();
130+
}
131+
rank = std::min((int)svd.singularValues().size(), rank);
132+
this->U_ = X.leftCols(n_cols_X)*svd.matrixU().leftCols(rank);
133+
this->Lambda_ = (svd.singularValues().head(rank).array().pow(2)-shift).matrix();
134+
return;
135+
}
136+
};
137+
138+
template<typename MatrixType>
139+
class REVD{
140+
private:
141+
std::unique_ptr<REVDStrategy<MatrixType>> revd_strategy_;
142+
DMatrix<double> U_;
143+
DVector<double> Lambda_;
144+
public:
145+
explicit REVD(std::unique_ptr<REVDStrategy<MatrixType>> &&strategy=std::make_unique<NysRSI<DMatrix<double>>>()): revd_strategy_(std::move(strategy)){}
146+
void compute(const MatrixType &A, int tr_rank, int max_iter=1e3){
147+
revd_strategy_->compute(A,tr_rank,max_iter);
148+
return;
149+
}
150+
//setters
151+
void setTol(double tol){ revd_strategy_->setTol(tol);}
152+
void setSeed(unsigned int seed){ revd_strategy_->setSeed(seed);}
153+
//getters
154+
int rank() const{ return revd_strategy_->rank();}
155+
DMatrix<double> matrixU() const{ return revd_strategy_->matrixU();}
156+
DVector<double> eigenValues() const{ return revd_strategy_->eigenValues();}
157+
};
158+
159+
}//core
160+
}//fdpade
161+
162+
#endif //REVD_H

0 commit comments

Comments
 (0)