The package originally started as just a header library to implement the sparse, symmetric matrix LDL decomposition algorithm described by Timothy Davis (https://fossies.org/linux/SuiteSparse/LDL/Doc/ldl_userguide.pdf), but now also includes a growing sparse matrix class that is compataible with the Eigen library. Sparse matrices are possible in Eigen, but hard to work with and lacks LDL decomposition, so this package is meant to provide a simple interface with relatively efficient functionality. However, these functions will be optimised in future versions. We use the sparse functionality in the glmmrBase and associated packages to store and use, for example, sparse design and covariance matrices.
The header file in /inst/include/sparse.h imports the classes for the library and can be included in other projects using Rcpp in R.
The R function sparse_chol() provides an R interface using compressed column form as per the CsparseMatrix in the Matrix package.
The sparse matrix class sparse can be initialised in several ways including from an Rcpp::NumericMatrix and an Eigen::MatrixXd. It is effectively a container for
a sparse matrix representation of the matrix. It can also just be initialised by specifying the rows and columns and later inserting elements. For example,
#include <sparse.h>
#include <RcppEigen.h> //assuming this is for an R package
// to create a 4x3 matrix, stored in row major order
sparse A(4,3,true);
A.insert(0,0,1);
A.insert(0,2,2);
A.insert(1,1,1);
A.insert(2,1,3);
A.insert(3,0,2);
A.insert(3,2,3);
// generate an Eigen matrix
MatrixXd B(3,3);
for(int i = 0; i < 3; i++){
for(int j = 0; j < 3; j++){
B(i,j) = i + 1 + j*3;
}
}
MatrixXd AB = A * B; // multiply sparse by sparse
// generate a transpose of the matrix
sparse At = A;
At.transpose();
// and multiply
At *= A;
The SparseChol class is responsible for the sparse LDL decomposition. We initialise the class using a sparse matrix.
#include <sparse.h>
SparseChol sp(mat); // where mat is of class sparse
int d = sp.ldl_numeric(); // run decomposition, returns dimension of matrix if success
//to retrieve the Cholesky decomposion LD^0.5
sparse LD = sp.LD();
// or we can return L and D separately
// to solve the linear system Ax = v where v is a std::vector<double> we pass a reference to the first element of the vector
// the system is solved in place so that v will contain the result
spchol.ldl_lsolve(&v[0]); // Lx = v
spchol.ldl_d2solve(&v[0]); // Dx = v
We can use the LDL decomposition in R with the associated functions provided for demonstration.
> n <- 10
> Ap <- c(0, 1, 2, 3, 4, 6, 7, 9, 11, 15, ANZ)
> Ai <- c(1, 2, 3, 4, 2,5, 6, 5,7, 5,8, 1,5,8,9, 2,5,7,10)
> Ax = c(1.7, 1., 1.5, 1.1, .02,2.6, 1.2, .16,1.3, .09,1.6,
.13,.52,.11,1.4, .01,.53,.56,3.1)
> Matrix::sparseMatrix(i = Ai, p=Ap, x=Ax, symmetric = TRUE)
10 x 10 sparse Matrix of class "dsCMatrix"
[1,] 1.70 . . . . . . . 0.13 .
[2,] . 1.00 . . 0.02 . . . . 0.01
[3,] . . 1.5 . . . . . . .
[4,] . . . 1.1 . . . . . .
[5,] . 0.02 . . 2.60 . 0.16 0.09 0.52 0.53
[6,] . . . . . 1.2 . . . .
[7,] . . . . 0.16 . 1.30 . . 0.56
[8,] . . . . 0.09 . . 1.60 0.11 .
[9,] 0.13 . . . 0.52 . . 0.11 1.40 .
[10,] . 0.01 . . 0.53 . 0.56 . . 3.10
> out <-sparse_chol(n,Ap,Ai,Ax)
> sparse_L(out)
10 x 10 sparse Matrix of class "dtCMatrix"
[1,] 1.00000000 . . . . . . . . .
[2,] . 1.00 . . . . . . . .
[3,] . . 1 . . . . . . .
[4,] . . . 1 . . . . . .
[5,] . 0.02 . . 1.00000000 . . . . .
[6,] . . . . . 1 . . . .
[7,] . . . . 0.06154793 . 1.000000000 . . .
[8,] . . . . 0.03462071 . -0.004293535 1.00000000 . .
[9,] 0.07647059 . . . 0.20003077 . -0.024807089 0.05752527 1.00000000 .
[10,] . 0.01 . . 0.20380058 . 0.408782664 -0.01006831 -0.07185228 1
> sparse_D(out)
10 x 10 diagonal matrix of class "ddiMatrix"
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 1.7 . . . . . . . . .
[2,] . 1 . . . . . . . .
[3,] . . 1.5 . . . . . . .
[4,] . . . 1.1 . . . . . .
[5,] . . . . 2.5996 . . . . .
[6,] . . . . . 1.2 . . . .
[7,] . . . . . . 1.290152 . . .
[8,] . . . . . . . 1.59686 . .
[9,] . . . . . . . . 1.279965 .
[10,] . . . . . . . . . 2.769568
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent