1%%%%%%%%%%%%%%%%%%% 2% XLiFE++ is an extended library of finite elements written in C++ 3% Copyright (C) 2014 Lunéville, Eric; Kielbasiewicz, Nicolas; Lafranche, Yvon; Nguyen, Manh-Ha; Chambeyron, Colin 4% 5% This program is free software: you can redistribute it and/or modify 6% it under the terms of the GNU General Public License as published by 7% the Free Software Foundation, either version 3 of the License, or 8% (at your option) any later version. 9% This program is distributed in the hope that it will be useful, 10% but WITHOUT ANY WARRANTY; without even the implied warranty of 11% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 12% GNU General Public License for more details. 13% You should have received a copy of the GNU General Public License 14% along with this program. If not, see <http://www.gnu.org/licenses/>. 15%%%%%%%%%%%%%%%%%%% 16 17\section{Approximate matrix} 18 19In the context of hierarchical matrix, some sub-matrices may be "compressed" to save memory and time computation. \xlifepp provides a general class {\class ApproximateMatrix} to deal with real or complex approximate matrix that can be 20\begin{itemize} 21 \item low rank matrix represented as the product $UDV^*$ ({\class LowRankMatrix}) 22 \item sinus cardinal decomposition ({\class SCDMatrix}), not yet available 23 \item fast multipole representation ({\class FMMMatrix}), not yet available 24\end{itemize} 25 26\subsection{The {\classtitle ApproximateMatrix} class} 27 28The {\class ApproximateMatrix} class is a template abstract class (templated by the type of coefficients) that manages only the type of approximation, the name of the approximate matrix and defines all virtual functions that children has to implement: 29\vspace{.1cm} 30\begin{lstlisting} 31template <typename T> 32class ApproximateMatrix 33{public : 34MatrixApproximationType approximationType; //type of approximation 35string_t name; //optional name useful for doc 36virtual ~ApproximateMatrix() {}; //virtual destructor 37virtual ApproximateMatrix<T>* clone() const=0; //creation of a clone 38virtual number_t numberOfRows() const =0; //nb of rows 39virtual number_t numberOfCols() const =0; //nb of columns 40virtual number_t nbNonZero() const =0; //nb of non zeros coefficient 41virtual number_t rank() const =0; //rank of matrix 42virtual vector<T>& multMatrixVector(const vector<T>&, vector<T>&) const = 0; 43virtual vector<T>& multVectorMatrix(const vector<T>&, vector<T>&) const = 0; 44virtual real_t norm2() const =0; //Frobenius norm 45virtual real_t norminfty() const=0; //infinite norm 46virtual LargeMatrix<T> toLargeMatrix() const = 0; //convert to LargeMatrix 47virtual void print(std::ostream&) const =0; //print to stream 48}; 49\end{lstlisting} 50\vspace{.1cm} 51Up to now, the only children available is the templated {\class LowRankMatrix<T>} class. 52 53\subsection{The {\classtitle LowRankMatrix} class} 54 55A low rank matrix is a matrix of the form 56$$U\,D\,V^*$$ 57where $U$ is a $m\times r$ matrix, $V$ is a $n\times r$ matrix and $D$ is a $r\times r$ diagonal matrix. The rank of a such matrix is at most $r$. So it is a low rank representation of a $m\times n$ matrix if $r$ is small compared to $m$, $n$. Contrary to the class name suggests, some matrix may be not low rank representation. Think about the SVD (singular value decomposition) of a $m\times n$ matrix that can be handled by this class ($r=\min(m,n)$).\\ 58 59The {\class LowRankMatrix} class is derived from the {\class ApproximateMatrix} and handles the $U,V$ matrices as {\class Matrix} objects of \xlifepp (dense matrix with major row access) and the diagonal matrix $D$ as a {\class Vector} object of \xlifepp : 60\vspace{.1cm} 61\begin{lstlisting} 62template <typename T> 63class LowRankMatrix : public ApproximateMatrix<T> 64{public : 65 Matrix<T> U_, V_; 66 Vector<T> D_; 67 ... 68} 69\end{lstlisting} 70\vspace{.1cm} 71 72Different constructors are available: 73\vspace{.1cm} 74\begin{lstlisting} 75LowRankMatrix(); //default constructor 76//constructor from explicit matrices and vector 77LowRankMatrix(const Matrix<T>&, const Matrix<T>&, const Vector<T>&, 78 const string_t& na=""); 79LowRankMatrix(const Matrix<T>&, const Matrix<T>&, const string_t& na=""); 80//constructor from dimensions 81LowRankMatrix(number_t m, number_t n, number_t r, const string_t& na=""); 82LowRankMatrix(number_t m, number_t n, number_t r, bool noDiag, 83 const string_t& na=""); 84//constructor from matrix iterators and vector iterator 85template <typename ITM, typename ITV> 86LowRankMatrix(number_t m, number_t n, number_t r, 87 ITM itmu, ITM itmv, ITV itd, const string_t& na=""); 88template <typename ITM> 89LowRankMatrix(number_t m, number_t n, number_t r, 90 ITM itmu, ITM itmv, const string_t& na=""); 91\end{lstlisting} 92\vspace{.1cm} 93The {\class LowRankMatrix} class provides all virtual functions declared by the {\class ApproximateMatrix} class: 94\vspace{.1cm} 95\begin{lstlisting} 96ApproximateMatrix<T>* clone() const; 97number_t numberOfRows() const; 98number_t numberOfCols() const; 99number_t nbNonZero() const; 100number_t rank() const; 101bool hasDiag() const; 102vector<T>& multMatrixVector(const vector<T>&, vector<T>&) const; 103vector<T>& multVectorMatrix(const vector<T>&, vector<T>&) const; 104LargeMatrix<T> toLargeMatrix() const; 105real_t norm2() const; 106real_t norminfty() const; 107void print(std::ostream&) const; 108\end{lstlisting} 109\vspace{.1cm} 110In addition, the class offers specific algebraic operations as member functions: 111\vspace{.1cm} 112\begin{lstlisting} 113LowRankMatrix<T> svd(real_t eps) const; 114LowRankMatrix<T>& operator+=(const LowRankMatrix<T>&); 115LowRankMatrix<T>& operator-=(const LowRankMatrix<T>&); 116LowRankMatrix<T>& operator*=(const T&); 117LowRankMatrix<T>& operator/=(const T&); 118LowRankMatrix<T>& add(const LowRankMatrix<T>&, const T&, real_t eps=0); 119\end{lstlisting} 120\vspace{.1cm} 121and as extern functions to the class (template declaration is omitted): 122\vspace{.1cm} 123\begin{lstlisting} 124LowRankMatrix<T> combine(const LowRankMatrix<T>&, const T&, 125 const LowRankMatrix<T>&, const T&, real_t eps=0); 126LowRankMatrix<T> operator+(const LowRankMatrix<T>&, const LowRankMatrix<T>&); 127LowRankMatrix<T> operator-(const LowRankMatrix<T>&, const LowRankMatrix<T>&); 128LowRankMatrix<T> operator-(const LowRankMatrix<T>&); 129LowRankMatrix<T> operator+(const LowRankMatrix<T>&); 130LowRankMatrix<T> operator*(const LowRankMatrix<T>&, const T&); 131LowRankMatrix<T> operator/(const LowRankMatrix<T>&, const T&); 132LowRankMatrix<T> operator*(const T&, const LowRankMatrix<T>&); 133\end{lstlisting} 134\vspace{.1cm} 135Note that the combination of a low rank matrix $L_1=U_1D_1V_1$ of rank $r_1$ and a low rank matrix $L_2=U_2D_2V_2$ of rank $r_2$ produces a low rank matrix $L=U\,D\,V$ of rank $r=r_1+r_2$ because 136$$L=a_1L_1+a_2L_2=\left[\begin{array}{cc}U_1 & U_2\end{array}\right] 137\left[\begin{array}{cc} 138a_1D_1 & 0 \\ 1390 & a_2D_2 140\end{array}\right] 141\left[\begin{array}{c}V_1\\ V_2\end{array}\right]. 142$$ 143When non zero, the {\tt eps} argument in the {\tt LowRankMatrix::add} and {\tt combine} functions is used to re-compress the combined matrix with a truncated svd.\\ 144 145Finally, some compression methods based on SVD are implemented too. These compression methods "produce" a {\class LowRankMatrix} from a dense matrix stored in major column access and passed by a pointer to its first coefficient (template declaration is omitted): 146\vspace{.1cm} 147\begin{lstlisting} 148void svdCompression(T* mat, number_t m, number_t n, number_t r, 149 LowRankMatrix<T>& lrm); 150void svdCompression(T* mat, number_t m, number_t n, real_t eps, 151 LowRankMatrix<T>& lrm); 152void rsvdCompression(T* mat, number_t m, number_t n, number_t r, 153 LowRankMatrix<T>& lrm); 154void rsvdCompression(T* mat, number_t m, number_t n, real_t eps, 155 LowRankMatrix<T>& lrm); 156void r3svdCompression(T* mat, number_t m, number_t n, real_t eps, 157 LowRankMatrix<T>& lrm, 158 number_t t = 0, number_t p = 0, number_t q = 0, 159 number_t maxit = 0) 160\end{lstlisting} 161\vspace{.1cm} 162The SVD compression functions do the full svd of the given matrix and then truncate the number of singular values and singular vectors keeping either the $r$ largest singular values/vectors or the singular values greater than a given {\tt eps}. The fact that this process leads to a good approximate matrix results from the Eckart–Young–Mirsky theorem.\\ 163 164As the full svd is an expansive algorithm time computation, some alternative methods based on the random svd are also available. Random svd consists in capturing the matrix range using only few gaussian random vectors and doing a svd on a smaller matrix. Again, two versions of random svd are available, one providing a low rank matrix of a given rank, the other one providing a low rank matrix with a control on the approximate error. This last method iterates on the rank, so it is more expansive. The r3svd is a more sophisticated iterative method (see https://arxiv.org/ftp/arxiv/papers/1605/1605.08134.pdf).\\ 165 166\begin{remark} 167As svd and qr algorithms have not been implemented for the {\class Matrix} class, \xlifepp uses the Eigen library that it is provided when installing \xlifepp. 168\end{remark} 169 170\begin{remark2} 171svd, rsvd and aca compression does not work for matrix of matrices! 172\end{remark2} 173 174\displayInfos{library=hierarchicalMatrix, header=ApproximateMatrix.hpp, implementation=ApproximateMatrix.cpp, test={test\_HMatrix.cpp}, header dep={largeMatrix.h, config.h, utils.h}} 175