Skylark (Sketching Library)  0.1
/var/lib/jenkins/jobs/Skylark/workspace/algorithms/regression/accelerated_linearl2_regression_solver_Elemental.hpp
Go to the documentation of this file.
00001 #ifndef SKYLARK_ACCELERATED_LINEARL2_REGRESSION_SOLVER_ELEMENTAL_HPP
00002 #define SKYLARK_ACCELERATED_LINEARL2_REGRESSION_SOLVER_ELEMENTAL_HPP
00003 
00004 #include <elemental.hpp>
00005 
00006 #include "regression_problem.hpp"
00007 
00008 namespace skylark {
00009 namespace algorithms {
00010 
00011 namespace flinl2_internal {
00012 
00013 extern "C" {
00014 
00015 void BLAS(dtrcon)(const char* norm, const char* uplo, const char* diag,
00016     const int* n, const double* A, const int *lda, double *rcond,
00017     double *work, int *iwork, int *info);
00018 
00019 void BLAS(strcon)(const char* norm, const char* uplo, const char* diag,
00020     const int* n, const float* A, const int *lda, float *rcond,
00021     float *work, int *iwork, int *info);
00022 
00023 }
00024 
00025 inline double utcondest(const elem::Matrix<double>& A) {
00026     int n = A.Height(), ld = A.LDim(), info;
00027     double *work = new double[3 * n];
00028     int *iwork = new int[n];
00029     double rcond;
00030     BLAS(dtrcon)("1", "U", "N", &n, A.LockedBuffer(), &ld, &rcond, work,
00031         iwork, &info); // TODO check info
00032     delete[] work;
00033     delete[] iwork;
00034     return 1/rcond;
00035 }
00036 
00037 inline double utcondest(const elem::Matrix<float>& A) {
00038     int n = A.Height(), ld = A.LDim(), info;
00039     float *work  = new float[3 * n];
00040     int *iwork = new int[n];
00041     float rcond;
00042     BLAS(strcon)("1", "U", "N", &n, A.LockedBuffer(), &ld, &rcond, work,
00043         iwork, &info);  // TODO check info
00044     delete[] work;
00045     delete[] iwork;
00046     return 1/rcond;
00047 }
00048 
00049 template<typename T>
00050 inline double utcondest(const elem::DistMatrix<T, elem::STAR, elem::STAR>& A) {
00051     return utcondest(A.LockedMatrix());
00052 }
00053 
00054 template<typename SolType, typename SketchType, typename PrecondType>
00055 double build_precond(SketchType& SA,
00056     PrecondType& R, nla::precond_t<SolType> *&P, qr_precond_tag) {
00057     elem::qr::Explicit(SA.Matrix(), R.Matrix()); // TODO
00058     P =
00059         new nla::tri_inverse_precond_t<SolType, PrecondType,
00060                                        elem::UPPER, elem::NON_UNIT>(R);
00061     return utcondest(R);
00062 }
00063 
00064 template<typename SolType, typename SketchType, typename PrecondType>
00065 double build_precond(SketchType& SA,
00066     PrecondType& V, nla::precond_t<SolType> *&P, svd_precond_tag) {
00067 
00068     int n = SA.Width();
00069     PrecondType s(SA);
00070     s.Resize(n, 1);
00071     elem::SVD(SA.Matrix(), s.Matrix(), V.Matrix()); // TODO
00072     for(int i = 0; i < n; i++)
00073         s.Set(i, 0, 1 / s.Get(i, 0));
00074     base::DiagonalScale(elem::RIGHT, elem::NORMAL, s, V);
00075     P =
00076         new nla::mat_precond_t<SolType, PrecondType>(V);
00077     return s.Get(0,0) / s.Get(n-1, 0);
00078 }
00079 
00080 }  // namespace flinl2_internal
00081 
00083 template <typename ValueType, elem::Distribution VD,
00084           template <typename, typename> class TransformType,
00085           typename PrecondTag>
00086 class accelerated_regression_solver_t<
00087     regression_problem_t<elem::DistMatrix<ValueType, VD, elem::STAR>,
00088                          linear_tag, l2_tag, no_reg_tag>,
00089     elem::DistMatrix<ValueType, VD, elem::STAR>,
00090     elem::DistMatrix<ValueType, elem::STAR, elem::STAR>,
00091     simplified_blendenpik_tag<TransformType, PrecondTag> > {
00092 
00093 public:
00094 
00095     typedef ValueType value_type;
00096 
00097     typedef elem::DistMatrix<ValueType, VD, elem::STAR> matrix_type;
00098     typedef elem::DistMatrix<ValueType, VD, elem::STAR> rhs_type;
00099     typedef elem::DistMatrix<ValueType, elem::STAR, elem::STAR> sol_type;
00100 
00101     typedef regression_problem_t<matrix_type,
00102                                  linear_tag, l2_tag, no_reg_tag> problem_type;
00103 
00104 private:
00105 
00106     typedef elem::DistMatrix<ValueType, elem::STAR, elem::STAR> precond_type;
00107     typedef precond_type sketch_type;
00108     // The assumption is that the sketch is not much bigger than the
00109     // preconditioner, so we should use the same matrix distribution.
00110 
00111     const int _m;
00112     const int _n;
00113     const matrix_type &_A;
00114     precond_type _R;
00115     nla::precond_t<sol_type> *_precond_R;
00116 
00117 public:
00123     accelerated_regression_solver_t(const problem_type& problem, base::context_t& context) :
00124         _m(problem.m), _n(problem.n), _A(problem.input_matrix),
00125         _R(_n, _n, problem.input_matrix.Grid()) {
00126         // TODO n < m ???
00127 
00128         int t = 4 * _n;    // TODO parameter.
00129 
00130         TransformType<matrix_type, sketch_type> S(_m, t, context);
00131         sketch_type SA(t, _n);
00132         S.apply(_A, SA, sketch::columnwise_tag());
00133 
00134         flinl2_internal::build_precond(SA, _R, _precond_R, PrecondTag());
00135     }
00136 
00137     ~accelerated_regression_solver_t() {
00138         delete _precond_R;
00139     }
00140 
00141     int solve(const rhs_type& b, sol_type& x) {
00142         return LSQR(_A, b, x, nla::iter_params_t(), *_precond_R);
00143     }
00144 };
00145 
00147 template <typename ValueType, elem::Distribution VD,
00148           typename PrecondTag>
00149 class accelerated_regression_solver_t<
00150     regression_problem_t<elem::DistMatrix<ValueType, VD, elem::STAR>,
00151                          linear_tag, l2_tag, no_reg_tag>,
00152     elem::DistMatrix<ValueType, VD, elem::STAR>,
00153     elem::DistMatrix<ValueType, elem::STAR, elem::STAR>,
00154     blendenpik_tag<PrecondTag> > {
00155 
00156 public:
00157 
00158     typedef ValueType value_type;
00159 
00160     typedef elem::DistMatrix<ValueType, VD, elem::STAR> matrix_type;
00161     typedef elem::DistMatrix<ValueType, VD, elem::STAR> rhs_type;
00162     typedef elem::DistMatrix<ValueType, elem::STAR, elem::STAR> sol_type;
00163 
00164     typedef regression_problem_t<matrix_type,
00165                                  linear_tag, l2_tag, no_reg_tag> problem_type;
00166 
00167 private:
00168 
00169     typedef elem::DistMatrix<ValueType, elem::STAR, elem::STAR> precond_type;
00170     typedef precond_type sketch_type; 
00171     // The assumption is that the sketch is not much bigger than the
00172     // preconditioner, so we should use the same matrix distribution.
00173 
00174     const int _m;
00175     const int _n;
00176     const matrix_type &_A;
00177     precond_type _R;
00178     nla::precond_t<sol_type> *_precond_R;
00179 
00180     regression_solver_t<problem_type, rhs_type, sol_type, svd_l2_solver_tag> 
00181     *_alt_solver;
00182 
00183 public:
00189     accelerated_regression_solver_t(const problem_type& problem, base::context_t& context) :
00190         _m(problem.m), _n(problem.n), _A(problem.input_matrix),
00191         _R(_n, _n, problem.input_matrix.Grid()) {
00192         // TODO n < m ???
00193 
00194         int t = 4 * _n;    // TODO parameter.
00195         double scale = std::sqrt((double)_m / (double)t);
00196 
00197         elem::DistMatrix<ValueType, elem::STAR, VD> Ar(_A.Grid());
00198         elem::DistMatrix<ValueType, elem::STAR, VD> dist_SA(t, _n, _A.Grid());
00199         sketch_type SA(t, _n, _A.Grid());
00200         boost::random::uniform_int_distribution<int> distribution(0, _m- 1);
00201 
00202         Ar = _A;
00203         double condest = 0;
00204         int attempts = 0;
00205         do {
00206             sketch::RFUT_t<elem::DistMatrix<ValueType, elem::STAR, VD>,
00207                            sketch::fft_futs<double>::DCT_t,
00208                            utility::rademacher_distribution_t<value_type> > 
00209                 F(_m, context);
00210             F.apply(Ar, Ar, sketch::columnwise_tag());
00211 
00212             std::vector<int> samples =
00213                 context.generate_random_samples_array(t, distribution);
00214             for (int j = 0; j < Ar.LocalWidth(); j++)
00215                 for (int i = 0; i < t; i++) {
00216                     int row = samples[i];
00217                     dist_SA.Matrix().Set(i, j, scale * Ar.Matrix().Get(row, j));
00218                 }
00219 
00220             SA = dist_SA;
00221             condest = flinl2_internal::build_precond(SA, _R, _precond_R, PrecondTag());
00222             attempts++;
00223         } while (condest > 1e14 && attempts < 3); // TODO parameters
00224 
00225         if (condest <= 1e14)
00226             _alt_solver = nullptr;
00227         else {
00228             _alt_solver =
00229                 new regression_solver_t<problem_type,
00230                                       rhs_type,
00231                                       sol_type, svd_l2_solver_tag>(problem);
00232             delete _precond_R;
00233             _precond_R = nullptr;
00234         }
00235     }
00236 
00237     ~accelerated_regression_solver_t() {
00238         if (_precond_R != nullptr)
00239             delete _precond_R;
00240         if (_alt_solver != nullptr)
00241             delete _alt_solver;
00242     }
00243 
00244     int solve(const rhs_type& b, sol_type& x) {
00245         return LSQR(_A, b, x, nla::iter_params_t(), *_precond_R);
00246     }
00247 };
00248 
00250 template <typename ValueType, elem::Distribution VD,
00251           typename PrecondTag>
00252 class accelerated_regression_solver_t<
00253     regression_problem_t<elem::DistMatrix<ValueType, VD, elem::STAR>,
00254                          linear_tag, l2_tag, no_reg_tag>,
00255     elem::DistMatrix<ValueType, VD, elem::STAR>,
00256     elem::DistMatrix<ValueType, elem::STAR, elem::STAR>,
00257     lsrn_tag<PrecondTag> > {
00258 
00259 public:
00260 
00261     typedef ValueType value_type;
00262 
00263     typedef elem::DistMatrix<ValueType, VD, elem::STAR> matrix_type;
00264     typedef elem::DistMatrix<ValueType, VD, elem::STAR> rhs_type;
00265     typedef elem::DistMatrix<ValueType, elem::STAR, elem::STAR> sol_type;
00266 
00267     typedef regression_problem_t<matrix_type,
00268                                  linear_tag, l2_tag, no_reg_tag> problem_type;
00269 
00270 private:
00271 
00272     typedef elem::DistMatrix<ValueType, elem::STAR, elem::STAR> precond_type;
00273     typedef precond_type sketch_type; 
00274     // The assumption is that the sketch is not much bigger than the
00275     // preconditioner, so we should use the same matrix distribution.
00276 
00277     const int _m;
00278     const int _n;
00279     const matrix_type &_A;
00280     bool _use_lsqr;
00281     double _sigma_U, _sigma_L;
00282     precond_type _R;
00283     nla::precond_t<sol_type> *_precond_R;
00284     nla::iter_params_t _params;
00285 
00286 public:
00292     accelerated_regression_solver_t(const problem_type& problem, base::context_t& context) :
00293         _m(problem.m), _n(problem.n), _A(problem.input_matrix),
00294         _R(_n, _n, problem.input_matrix.Grid()) {
00295         // TODO n < m ???
00296 
00297         int t = 4 * _n;    // TODO parameter.
00298         double epsilon = 1e-14;  // TODO parameter
00299         double delta = 1e-6; // TODO parameter
00300 
00301         sketch::JLT_t<matrix_type, sketch_type> S(_m, t, context);
00302         sketch_type SA(t, _n);
00303         S.apply(_A, SA, sketch::columnwise_tag());
00304         flinl2_internal::build_precond(SA, _R, _precond_R, PrecondTag());
00305 
00306         // Select alpha so that probability of failure is delta.
00307         // If alpha is too big, we need to use LSQR (although ill-conditioning
00308         // might be very severe so to prevent convergence).
00309         double alpha = std::sqrt(2 * std::log(2.0 / delta) / t);
00310         if (alpha >= (1 - std::sqrt(_n / t)))
00311             _use_lsqr = true;
00312         else {
00313             _use_lsqr = false;
00314             _sigma_U = std::sqrt(t) / ((1 - alpha) * std::sqrt(t)
00315                 - std::sqrt(_n));
00316             _sigma_L = std::sqrt(t) / ((1 + alpha) * std::sqrt(t)
00317                 + std::sqrt(_n));
00318         }
00319     }
00320 
00321     ~accelerated_regression_solver_t() {
00322         delete _precond_R;
00323     }
00324 
00325     int solve(const rhs_type& b, sol_type& x) {
00326         int ret;
00327         if (_use_lsqr)
00328             ret = LSQR(_A, b, x, _params, *_precond_R);
00329         else {
00330             ChebyshevLS(_A, b, x,  _sigma_L, _sigma_U,
00331                 _params, *_precond_R);
00332             ret = -6; // TODO! - check!
00333         }
00334         return ret; // TODO!
00335     }
00336 };
00337 
00338 
00339 
00340 } } 
00342 #endif // SKYLARK_ACCELERATED_LINEARL2_REGRESSION_SOLVER_ELEMENTAL_HPP