From b30c440ee789d2b8746ef5d65a68a0dba11df40d Mon Sep 17 00:00:00 2001 From: Yuriel Date: Fri, 9 May 2025 16:06:34 +0200 Subject: [PATCH 1/8] implement fast rrLU --- include/xfac/matrix/mat_decomp.h | 36 ++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/include/xfac/matrix/mat_decomp.h b/include/xfac/matrix/mat_decomp.h index 4c95f6b..c9da6e0 100644 --- a/include/xfac/matrix/mat_decomp.h +++ b/include/xfac/matrix/mat_decomp.h @@ -148,6 +148,42 @@ struct RRLUDecomp { readLU(A); } + void calculateFast(arma::Mat A, double reltol, int rankMax) + { + npivot= std::min(A.n_rows, A.n_cols); + if (reltol==0) reltol=npivot*std::abs(std::numeric_limits::epsilon()); + if (rankMax>0 && rankMax Q,R; + arma::qr(Q, R, J, A, "vector"); + arma::vec sv=R.diag(); + npivot=arma::find(sv>reltol*sv[0]).eval().size(); + J=J.head_cols(npivot); + } + arma::uvec I(npivot); + { // row permutation by LU + arma::Mat P; // P*A==L*U + lu(L, U, P, A.cols(J)); + for(auto i=0u; i D=L.diag().eval(); + L *= arma::diagmat(arma::inv(D)); + U *= arma::diagmat(D); + } + Iset=arma::conv_to>::from(I); + Jset=arma::conv_to>::from(J); + } + vector row_pivots() const { return {Iset.begin(), Iset.begin()+npivot}; } vector col_pivots() const { return {Jset.begin(), Jset.begin()+npivot}; } arma::Mat PivotMatrixTri() const From 2d61c9f2a6e2fb569d1cdfc443caa56202844244 Mon Sep 17 00:00:00 2001 From: Yuriel Date: Wed, 14 May 2025 10:15:01 +0200 Subject: [PATCH 2/8] fix previous --- include/xfac/matrix/mat_decomp.h | 106 ++++++++++++++++--------------- 1 file changed, 56 insertions(+), 50 deletions(-) diff --git a/include/xfac/matrix/mat_decomp.h b/include/xfac/matrix/mat_decomp.h index c9da6e0..029f770 100644 --- a/include/xfac/matrix/mat_decomp.h +++ b/include/xfac/matrix/mat_decomp.h @@ -2,7 +2,6 @@ #define MAT_DECOMP_H #include "xfac/index_set.h" -#include #define ARMA_DONT_USE_OPENMP #include @@ -96,6 +95,37 @@ struct MatSVDFixedTol: public MatDecompFixedTol> { //--------------------------------------- rank-revealing LU decomposition ----------------------- + +/// Given the pivot matrix P already in LU form, update the C columns according to P. +template +void apply_LU_on_cols(arma::Mat& C, arma::Mat const& P, bool leftOrthogonal) +{ + if (P.n_rows != C.n_rows) + throw std::invalid_argument("apply_LU_on_Cols: incompatible matrices"); + for(auto k=0u; k +void apply_LU_on_rows(arma::Mat& R, arma::Mat const& P, bool leftOrthogonal) +{ + if (P.n_cols != R.n_cols) + throw std::invalid_argument("apply_LU_on_Rows: incompatible matrices"); + for(auto k=0u; k struct RRLUDecomp { using value_type=T; @@ -154,32 +184,37 @@ struct RRLUDecomp { if (reltol==0) reltol=npivot*std::abs(std::numeric_limits::epsilon()); if (rankMax>0 && rankMax Q,R; arma::qr(Q, R, J, A, "vector"); - arma::vec sv=R.diag(); - npivot=arma::find(sv>reltol*sv[0]).eval().size(); - J=J.head_cols(npivot); - } - arma::uvec I(npivot); - { // row permutation by LU - arma::Mat P; // P*A==L*U - lu(L, U, P, A.cols(J)); - for(auto i=0u; ireltol*sv[0]).eval().size(); + if (rank P; + lu(L, U, P, A.cols(J.head_rows(npivot))); + for(auto i=0u; i D=U.diag().eval(); if (!leftOrthogonal) { - arma::Col D=L.diag().eval(); - L *= arma::diagmat(arma::inv(D)); - U *= arma::diagmat(D); + L = L*arma::diagmat(D); + U = arma::diagmat(1.0/D)*U; + } + if (npivot P=L.head_rows(npivot)+U; // pivot matrix in LU form + P.diag()=D; + apply_LU_on_cols(U2,P,leftOrthogonal); + U=arma::join_horiz(U,U2).eval(); } + Iset=arma::conv_to>::from(I); Jset=arma::conv_to>::from(J); } @@ -227,35 +262,6 @@ struct RRLUDecomp { } }; -/// Given the pivot matrix P already in LU form, update the C columns according to P. -template -void apply_LU_on_cols(arma::Mat& C, arma::Mat const& P, bool leftOrthogonal) -{ - if (P.n_rows != C.n_rows) - throw std::invalid_argument("apply_LU_on_Cols: incompatible matrices"); - for(auto k=0u; k -void apply_LU_on_rows(arma::Mat& R, arma::Mat const& P, bool leftOrthogonal) -{ - if (P.n_cols != R.n_cols) - throw std::invalid_argument("apply_LU_on_Rows: incompatible matrices"); - for(auto k=0u; k struct MatRRLUFixedTol: public MatDecompFixedTol> { @@ -294,7 +300,7 @@ struct ARRLUDecomp: public RRLUDecomp { RRLUDecomp::Iset=iota(fA.n_rows); RRLUDecomp::Jset=iota(fA.n_cols); - if (param.fullPiv) { this->calculate(fA.submat(Iset,Jset), param.reltol, param.bondDim); return; } + if (param.fullPiv) { this->calculateFast(fA.submat(Iset,Jset), param.reltol, param.bondDim); return; } int rankMax=std::min(fA.n_rows,fA.n_cols); if (param.bondDim>0 && param.bondDim Date: Wed, 14 May 2025 10:15:15 +0200 Subject: [PATCH 3/8] add a test for armadillo --- test/test_main.cpp | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/test/test_main.cpp b/test/test_main.cpp index afeca98..f5725ae 100644 --- a/test/test_main.cpp +++ b/test/test_main.cpp @@ -1,2 +1,20 @@ #define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one cpp file #include +#include + + +TEST_CASE("Test arma") { + using namespace arma; + mat A(120,10000,fill::randn), L, U, P; + lu(L, U, P, A); + + REQUIRE(norm(P*A-L*U)<1e-14*norm(A)); + std::cout << norm(L.diag()-1.0) << std::endl; + + auto D=U.diag().eval(); + L = L*arma::diagmat(D); + U = arma::diagmat(1.0/D)*U; + + REQUIRE(norm(P*A-L*U)<1e-14*norm(A)); + std::cout << norm(U.diag()-1.0) << std::endl; +} From 73df7b01f6fbdf42981547bed82b67dc300f1cc6 Mon Sep 17 00:00:00 2001 From: Yuriel Date: Thu, 15 May 2025 11:59:02 +0200 Subject: [PATCH 4/8] add the option to use the fast rrlu --- include/xfac/matrix/mat_decomp.h | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/include/xfac/matrix/mat_decomp.h b/include/xfac/matrix/mat_decomp.h index 029f770..572a87f 100644 --- a/include/xfac/matrix/mat_decomp.h +++ b/include/xfac/matrix/mat_decomp.h @@ -138,11 +138,11 @@ struct RRLUDecomp { RRLUDecomp()=default; - RRLUDecomp(arma::Mat const& A, bool leftOrthogonal_=true, double reltol=0, int rankMax=0) + RRLUDecomp(arma::Mat const& A, bool leftOrthogonal_=true, double reltol=0, int rankMax=0, bool fast=false) : leftOrthogonal(leftOrthogonal_) , Iset (iota(A.n_rows)) , Jset (iota(A.n_cols)) - { calculate(A, reltol, rankMax); } + { if (fast) calculateFast(A,reltol,rankMax); else calculate(A, reltol, rankMax); } void calculate(arma::Mat A, double reltol, int rankMax) { @@ -286,7 +286,7 @@ struct ARRLUParam { }; template -struct ARRLUDecomp: public RRLUDecomp { +struct ARRLUDecomp: public RRLUDecomp { //TODO: try to use the fast rrlu using RRLUDecomp::Iset; using RRLUDecomp::Jset; using RRLUDecomp::npivot; @@ -362,8 +362,8 @@ struct CURDecomp: public ARRLUDecomp { const arma::Mat C, R; - CURDecomp(arma::Mat const& A_, bool leftOrthogonal_=true, double reltol=1e-14, int rankMax=0) - : ARRLUDecomp(A_,leftOrthogonal_,reltol,rankMax) + CURDecomp(arma::Mat const& A_, bool leftOrthogonal_=true, double reltol=1e-14, int rankMax=0, bool fast=false) + : ARRLUDecomp(A_,leftOrthogonal_,reltol,rankMax,fast) , C(A_.cols(arma::conv_to::from(col_pivots()))) , R(A_.rows(arma::conv_to::from(row_pivots()))) {} From 83e6ddee51e4176c7584433b987cfe34d2de2f0b Mon Sep 17 00:00:00 2001 From: Yuriel Date: Thu, 15 May 2025 11:59:33 +0200 Subject: [PATCH 5/8] test the timing for the new rrlu --- test/test_matrix_ci.cpp | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/test/test_matrix_ci.cpp b/test/test_matrix_ci.cpp index 10f381a..23c5367 100644 --- a/test/test_matrix_ci.cpp +++ b/test/test_matrix_ci.cpp @@ -2,6 +2,7 @@ #include "xfac/tensor/tensor_ci.h" #include "xfac/grid.h" #include "xfac/matrix/mat_decomp.h" +#include using namespace arma; using namespace xfac; @@ -179,6 +180,33 @@ TEST_CASE("cur") CURDecomp sol(A); REQUIRE(arma::norm(A-sol.left()*sol.right()) sol(A,true,0.0,0); + auto t1 = std::chrono::high_resolution_clock::now(); + auto dt=std::chrono::duration_cast(t1-t0).count(); + cout<<" "< sol(A,true,0.0,0,fast); + auto t1 = std::chrono::high_resolution_clock::now(); + auto dt=std::chrono::duration_cast(t1-t0).count(); + cout<<" "< Date: Thu, 28 Aug 2025 10:19:09 +0200 Subject: [PATCH 6/8] add python binding to matrix decomposition by CI --- python/python_bindings.cpp | 13 ++++++++++++- 1 file changed, 12 insertions(+), 1 deletion(-) diff --git a/python/python_bindings.cpp b/python/python_bindings.cpp index 5a8b6ee..4487240 100644 --- a/python/python_bindings.cpp +++ b/python/python_bindings.cpp @@ -24,11 +24,22 @@ using cmpx=complex; template void declare_TensorCI(py::module &m, std::string typestr) { + + m.def(("decomposeCI"s+typestr).c_str(),[](arma::Mat A,bool isleft,double reltol,int maxBondDim) + { + // A.print("A"); + arma::Mat a=A;//carma::arr_to_mat(A); + // return; + auto ci=MatCURFixedTol {reltol, maxBondDim}; + auto [L,R]=ci(a,isleft); + return std::make_pair(L,R); + }, "A"_a, "leftCanonic"_a=true, "reltol"_a=1e-12, "maxBondDim"_a=0); + using TensorTraind=TensorTrain; py::class_(m, ("TensorTrain"s+typestr).c_str()) .def(py::init(), "len"_a) .def_readonly("core", &TensorTraind::M) - .def("setCoreAt",[](TensorTraind &tt, int i, py::array const& Mi){ tt.M.at(i)=carma::arr_to_cube(Mi); }, "i"_a, "Mi"_a) + .def("setCoreAt",[](TensorTraind &tt, int i, arma::Cube Mi){ tt.M.at(i)=Mi; }, "i"_a, "Mi"_a) .def("eval", &TensorTraind::eval) .def("sum", &TensorTraind::sum) .def("sum1", &TensorTraind::sum1) From 92cac26bf773282087f68ff5c4c04351a9d6a42e Mon Sep 17 00:00:00 2001 From: Yuriel Date: Thu, 28 Aug 2025 10:27:34 +0200 Subject: [PATCH 7/8] update readme --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index c80b528..b73d566 100644 --- a/README.md +++ b/README.md @@ -74,7 +74,7 @@ If you want: * To compile the **tests**, add `-D XFAC_BUILD_TEST=ON` to the cmake line above. -* To generate the **python** library, add `-D XFAC_BUILD_PYTHON=ON` +* To generate the **python** library, add `-D XFAC_BUILD_PYTHON=ON`. If you find any problem try with an older numpy `pip install "numpy<2.0"`. * To **install** the library in specific **local directory**, add `-D CMAKE_INSTALL_PREFIX=` From 8a98d90d67ef8fc4a92828e234a648db83a943a8 Mon Sep 17 00:00:00 2001 From: Yuriel Date: Fri, 29 Aug 2025 12:00:00 +0200 Subject: [PATCH 8/8] add argument to python decomposeCI --- python/python_bindings.cpp | 12 ++++-------- 1 file changed, 4 insertions(+), 8 deletions(-) diff --git a/python/python_bindings.cpp b/python/python_bindings.cpp index 4487240..8ecf1d5 100644 --- a/python/python_bindings.cpp +++ b/python/python_bindings.cpp @@ -25,15 +25,11 @@ using cmpx=complex; template void declare_TensorCI(py::module &m, std::string typestr) { - m.def(("decomposeCI"s+typestr).c_str(),[](arma::Mat A,bool isleft,double reltol,int maxBondDim) + m.def(("decomposeCI"s+typestr).c_str(),[](arma::Mat A,bool isleft,double reltol,int maxBondDim,bool fast) { - // A.print("A"); - arma::Mat a=A;//carma::arr_to_mat(A); - // return; - auto ci=MatCURFixedTol {reltol, maxBondDim}; - auto [L,R]=ci(a,isleft); - return std::make_pair(L,R); - }, "A"_a, "leftCanonic"_a=true, "reltol"_a=1e-12, "maxBondDim"_a=0); + auto ci=CURDecomp(A,isleft,reltol,maxBondDim,fast); + return std::make_pair(ci.left(),ci.right()); + }, "A"_a, "leftCanonic"_a=true, "reltol"_a=1e-12, "maxBondDim"_a=0, "fast"_a=true); using TensorTraind=TensorTrain; py::class_(m, ("TensorTrain"s+typestr).c_str())