48 #include "ROL_LinearAlgebra.hpp"
49 #include "ROL_LAPACK.hpp"
68 std::vector<Ptr<Vector<Real>>>
Y_;
88 Ptr<Elementwise::NormalRandom<Real>>
nrand_;
90 Ptr<std::normal_distribution<Real>>
dist_;
93 const int nvec(Y.size());
94 const Real
zero(0), one(1);
96 std::vector<Real> normQ(nvec,0);
98 for (
int j = 0; j < nvec; ++j) {
100 if (rjj > ROL_EPSILON<Real>()) {
101 for (
int k = 0; k <
orthIt_; ++k) {
102 for (
int i = 0; i < j; ++i) {
103 rij = Y[i]->dot(*Y[j]);
104 Y[j]->axpy(-rij,*Y[i]);
106 normQ[j] = Y[j]->norm();
108 for (
int i = 0; i < j; ++i) {
109 rij = std::abs(Y[i]->dot(*Y[j]));
110 if (rij >
orthTol_*normQ[j]*normQ[i]) {
119 if (rjj >
zero) Y[j]->scale(one/rjj);
123 int LSsolver(LA::Matrix<Real> &A, LA::Matrix<Real> &B,
bool trans =
false)
const {
125 char TRANS = (trans ?
'T' :
'N');
128 int NRHS = B.numCols();
130 int LDB = std::max(M,N);
131 std::vector<Real> WORK(1);
134 lapack_.GELS(TRANS,M,N,NRHS,A.values(),LDA,B.values(),LDB,&WORK[0],LWORK,&INFO);
136 LWORK =
static_cast<int>(WORK[0]);
138 lapack_.GELS(TRANS,M,N,NRHS,A.values(),LDA,B.values(),LDB,&WORK[0],LWORK,&INFO);
149 int K = std::min(M,N);
151 std::vector<Real> S(K);
152 LA::Matrix<Real> U(M,K);
154 LA::Matrix<Real> VT(K,N);
156 std::vector<Real> WORK(1), WORK0(1);
159 lapack_.GESVD(JOBU,JOBVT,M,N,A.values(),LDA,&S[0],U.values(),LDU,VT.values(),LDVT,&WORK[0],LWORK,&WORK0[0],&INFO);
160 LWORK =
static_cast<int>(WORK[0]);
162 lapack_.GESVD(JOBU,JOBVT,M,N,A.values(),LDA,&S[0],U.values(),LDU,VT.values(),LDVT,&WORK[0],LWORK,&WORK0[0],&INFO);
163 for (
int i = 0; i < M; ++i) {
164 for (
int j = 0; j < N; ++j) {
166 for (
int k = 0; k < r; ++k) {
167 A(i,j) += S[k] * U(i,k) * VT(k,j);
180 int K = std::min(M,N);
182 std::vector<Real> TAU(K);
183 std::vector<Real> WORK(1);
186 lapack_.GEQRF(M,N,
X_.values(),LDA,&TAU[0],&WORK[0],LWORK,&INFO);
187 LWORK =
static_cast<int>(WORK[0]);
189 lapack_.GEQRF(M,N,
X_.values(),LDA,&TAU[0],&WORK[0],LWORK,&INFO);
192 lapack_.ORGQR(M,N,K,
X_.values(),LDA,&TAU[0],&WORK[0],LWORK,&INFO);
193 LWORK =
static_cast<int>(WORK[0]);
195 lapack_.ORGQR(M,N,K,
X_.values(),LDA,&TAU[0],&WORK[0],LWORK,&INFO);
210 int infoP(0), infoQ(0), infoLS1(0), infoLS2(0), infoLRA(0);
215 LA::Matrix<Real> L(
s_,
k_), R(
s_,
k_);
216 for (
int i = 0; i <
s_; ++i) {
217 for (
int j = 0; j <
k_; ++j) {
218 L(i,j) =
Phi_[i]->dot(*
Y_[j]);
220 for (
int k = 0; k <
ncol_; ++k) R(i,j) +=
Psi_(k,i) *
X_(k,j);
225 LA::Matrix<Real> Zmat(s_,
k_);
226 for (
int i = 0; i <
k_; ++i) {
227 for (
int j = 0; j <
s_; ++j) Zmat(j,i) =
Z_(i,j);
230 for (
int i = 0; i <
k_; ++i) {
231 for (
int j = 0; j <
k_; ++j)
C_(j,i) = Zmat(i,j);
238 return std::abs(infoP)+std::abs(infoQ)+std::abs(infoLS1)
239 +std::abs(infoLS2)+std::abs(infoLRA);
246 Real orthTol = 1e-8,
int orthIt = 2,
bool truncate =
false,
247 unsigned dom_seed = 0,
unsigned rng_seed = 0)
252 nrand_ = makePtr<Elementwise::NormalRandom<Real>>(mu,sig,dom_seed);
253 unsigned seed = rng_seed;
254 if (seed == 0) seed = std::chrono::system_clock::now().time_since_epoch().count();
255 gen_ = makePtr<std::mt19937_64>(seed);
256 dist_ = makePtr<std::normal_distribution<Real>>(mu,sig);
259 rank_ = std::min(rank, maxRank_);
260 k_ = std::min(2*rank_+1, maxRank_);
261 s_ = std::min(2*
k_ +1, maxRank_);
265 for (
int i = 0; i <
k_; ++i) {
279 X_.putScalar(zero);
Z_.putScalar(zero);
C_.putScalar(zero);
280 for (
int i = 0; i <
k_; ++i)
Y_[i]->
zero();
283 for (
int i = 0; i <
s_; ++i) {
285 for (
int j = 0; j <
ncol_; ++j)
Psi_(j,i) = (*dist_)(*gen_);
287 for (
int i = 0; i <
k_; ++i) {
289 for (
int j = 0; j <
ncol_; ++j)
Omega_(j,i) = (*dist_)(*gen_);
297 int sold =
s_, kold =
k_;
303 for (
int i = sold; i <
s_; ++i)
Phi_.push_back(
Phi_[0]->clone());
306 for (
int i = kold; i <
k_; ++i) {
307 Y_.push_back(
Y_[0]->clone());
312 if (
out_ != nullPtr ) {
313 *
out_ << std::string(80,
'=') << std::endl;
314 *
out_ <<
" ROL::Sketch::setRank" << std::endl;
315 *
out_ <<
" **** Rank = " <<
rank_ << std::endl;
316 *
out_ <<
" **** k = " <<
k_ << std::endl;
317 *
out_ <<
" **** s = " <<
s_ << std::endl;
318 *
out_ << std::string(80,
'=') << std::endl;
328 if ( col >=
ncol_ || col < 0 )
return 1;
330 for (
int i = 0; i <
k_; ++i) {
332 for (
int j = 0; j <
ncol_; ++j)
X_(j,i) *= eta;
340 for (
int i = 0; i <
s_; ++i) {
342 for (
int j = 0; j <
s_; ++j) {
344 Z_(i,j) += nu*
Psi_(col,j)*hphi;
347 if (
out_ != nullPtr ) {
348 *
out_ << std::string(80,
'=') << std::endl;
349 *
out_ <<
" ROL::Sketch::advance" << std::endl;
350 *
out_ <<
" **** col = " << col << std::endl;
351 *
out_ <<
" **** norm(h) = " << h.
norm() << std::endl;
352 *
out_ << std::string(80,
'=') << std::endl;
364 if ( col >=
ncol_ || col < 0 )
return 2;
369 if (flag > 0 )
return 3;
372 if (flag > 0 )
return 4;
375 if (flag > 0 )
return 5;
379 for (
int i = 0; i <
k_; ++i) {
381 for (
int j = 0; j <
k_; ++j) coeff +=
C_(i,j) *
X_(col,j);
384 if (
out_ != nullPtr ) {
385 *
out_ << std::string(80,
'=') << std::endl;
386 *
out_ <<
" ROL::Sketch::reconstruct" << std::endl;
387 *
out_ <<
" **** col = " << col << std::endl;
388 *
out_ <<
" **** norm(a) = " << a.
norm() << std::endl;
389 *
out_ << std::string(80,
'=') << std::endl;
394 bool test(
const int rank, std::ostream &outStream = std::cout,
const int verbosity = 0) {
395 const Real one(1), tol(std::sqrt(ROL_EPSILON<Real>()));
396 using seed_type = std::mt19937_64::result_type;
397 seed_type
const seed = 123;
398 std::mt19937_64 eng{seed};
399 std::uniform_real_distribution<Real> dist(static_cast<Real>(0),static_cast<Real>(1));
401 std::vector<Ptr<Vector<Real>>> U(rank);
402 LA::Matrix<Real>
V(
ncol_,rank);
403 for (
int i = 0; i < rank; ++i) {
404 U[i] =
Y_[0]->clone();
405 U[i]->randomize(-one,one);
406 for (
int j = 0; j <
ncol_; ++j)
V(j,i) = dist(eng);
410 std::vector<Ptr<Vector<Real>>> A(
ncol_);
411 for (
int i = 0; i <
ncol_; ++i) {
412 A[i] =
Y_[0]->clone(); A[i]->zero();
413 for (
int j = 0; j < rank; ++j) {
414 A[i]->axpy(
V(i,j),*U[j]);
419 bool flagP =
testP(outStream, verbosity);
421 bool flagQ =
testQ(outStream, verbosity);
423 Real nerr(0), maxerr(0);
424 Ptr<Vector<Real>> err =
Y_[0]->clone();
425 for (
int i = 0; i <
ncol_; ++i) {
427 err->axpy(-one,*A[i]);
429 maxerr = (nerr > maxerr ? nerr : maxerr);
432 std::ios_base::fmtflags oflags(outStream.flags());
433 outStream << std::scientific << std::setprecision(3) << std::endl;
434 outStream <<
" TEST RECONSTRUCTION: Max Error = "
435 << std::setw(12) << std::right << maxerr
436 << std::endl << std::endl;
437 outStream.flags(oflags);
439 return flagP & flagQ & (maxerr < tol ?
true :
false);
445 bool testQ(std::ostream &outStream = std::cout,
const int verbosity = 0) {
446 const Real one(1), tol(std::sqrt(ROL_EPSILON<Real>()));
448 Real qij(0), err(0), maxerr(0);
449 std::ios_base::fmtflags oflags(outStream.flags());
450 if (verbosity > 0) outStream << std::scientific << std::setprecision(3);
452 outStream << std::endl
453 <<
" Printing Q'Q...This should be approximately equal to I"
454 << std::endl << std::endl;
456 for (
int i = 0; i <
k_; ++i) {
457 for (
int j = 0; j <
k_; ++j) {
458 qij =
Y_[i]->dot(*
Y_[j]);
459 err = (i==j ? std::abs(qij-one) : std::abs(qij));
460 maxerr = (err > maxerr ? err : maxerr);
461 if (verbosity > 1) outStream << std::setw(12) << std::right << qij;
463 if (verbosity > 1) outStream << std::endl;
464 if (maxerr > tol)
break;
467 outStream << std::endl <<
" TEST ORTHOGONALIZATION: Max Error = "
468 << std::setw(12) << std::right << maxerr
470 outStream.flags(oflags);
472 return (maxerr < tol ?
true :
false);
475 bool testP(std::ostream &outStream = std::cout,
const int verbosity = 0) {
476 const Real
zero(0), one(1), tol(std::sqrt(ROL_EPSILON<Real>()));
478 Real qij(0), err(0), maxerr(0);
479 std::ios_base::fmtflags oflags(outStream.flags());
480 if (verbosity > 0) outStream << std::scientific << std::setprecision(3);
482 outStream << std::endl
483 <<
" Printing P'P...This should be approximately equal to I"
484 << std::endl << std::endl;
486 for (
int i = 0; i <
k_; ++i) {
487 for (
int j = 0; j <
k_; ++j) {
489 for (
int k = 0; k <
ncol_; ++k) qij +=
X_(k,i) *
X_(k,j);
490 err = (i==j ? std::abs(qij-one) : std::abs(qij));
491 maxerr = (err > maxerr ? err : maxerr);
492 if (verbosity > 1) outStream << std::setw(12) << std::right << qij;
494 if (verbosity > 1) outStream << std::endl;
495 if (maxerr > tol)
break;
498 outStream << std::endl <<
" TEST ORTHOGONALIZATION: Max Error = "
499 << std::setw(12) << std::right << maxerr
501 outStream.flags(oflags);
503 return (maxerr < tol ?
true :
false);
void setStream(Ptr< std::ostream > &out)
Ptr< std::normal_distribution< Real > > dist_
std::vector< Ptr< Vector< Real > > > Y_
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
int reconstruct(Vector< Real > &a, const int col)
Ptr< std::mt19937_64 > gen_
Contains definitions of custom data types in ROL.
Provides an interface for randomized sketching.
int LSsolver(LA::Matrix< Real > &A, LA::Matrix< Real > &B, bool trans=false) const
std::vector< Ptr< Vector< Real > > > Upsilon_
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
Ptr< Elementwise::NormalRandom< Real > > nrand_
virtual Real dot(const Vector &x) const =0
Compute where .
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Sketch(const Vector< Real > &x, int ncol, int rank, Real orthTol=1e-8, int orthIt=2, bool truncate=false, unsigned dom_seed=0, unsigned rng_seed=0)
bool testQ(std::ostream &outStream=std::cout, const int verbosity=0)
void reset(bool randomize=true)
void mgs2(std::vector< Ptr< Vector< Real >>> &Y) const
bool testP(std::ostream &outStream=std::cout, const int verbosity=0)
bool test(const int rank, std::ostream &outStream=std::cout, const int verbosity=0)
LAPACK< int, Real > lapack_
LA::Matrix< Real > Omega_
int advance(Real nu, const Vector< Real > &h, int col, Real eta=1.0)
virtual Real norm() const =0
Returns where .
std::vector< Ptr< Vector< Real > > > Phi_
int lowRankApprox(LA::Matrix< Real > &A, int r) const