STK++ 0.9.13
Use of the Algebra utilities

Introduction to Algebra project

The Algebra project provides structures, tools and methods of the usual algebra techniques on matrices and vectors. If lapack is installed on your system, most of the tools presented here use it in background. Otherwise built-in (but maybe less efficient) methods are used.

Main operations that can be performed are:

Matrix decomposition

All methods for matrix decomposition are enclosed in classes. Computations are launched using the run method.

QR decomposition

QR decomposition (also called a QR factorization) of a matrix is a decomposition of a matrix A into a product A = QR of an orthogonal matrix Q and an upper triangular matrix R. Given a matrix A of size (m,n), QR decomposition (also called a QR factorization) of A is a decomposition into a product A=QR of an orthogonal matrix Q of size (m,m) and an upper triangular matrix R of size (m,n).

QR decomposition can be achieved using either the class

STK::Qr
STK::Qr
The class Qr perform the QR decomposition of an ArrayXX.
Definition STK_Qr.h:98

or if lapack is available

STK::lapack::Qr

In the later case, your code have to be compiled using -DSTKUSELAPACK flag and linked using -llapack (for GNU-like compilers).

ExampleOutput
#include <STKpp.h>
using namespace STK;
int main(int argc, char** argv)
{
Array2D<Real> a(5,4);
a << 0, 1, 2, 3,
2, 3, 4, 5,
2, 1, 6, 7,
0, 3,-1, 2,
3,-1, 1, 1;
stk_cout << _T("STK++ QR decomposition:\n");
stk_cout << _T("-----------------------\n");
Qr q(a); q.run();
stk_cout << _T("R matrix:\n");
stk_cout << q.R();
ArrayXX QR = q.R();
applyLeftHouseholderArray(QR, q.Q());
stk_cout << _T("QR matrix:\n");
stk_cout << QR;
stk_cout << _T("\nlapack QR decomposition:\n");
stk_cout << _T("--------------------------\n");
lapack::Qr lq(a); lq.run();
stk_cout << _T("R matrix:\n");
stk_cout << lq.R();
QR = lq.R();
applyLeftHouseholderArray(QR, lq.Q());
stk_cout << _T("QR matrix:\n");
stk_cout << QR;
return 0;
}
stk_cout
#define stk_cout
Standard stk output stream.
Definition STK_Typedefs.h:73
_T
#define _T(x)
Let x unmodified.
Definition STK_Typedefs.h:71
STKpp.h
This file include all the header files of the STK++ project.
STK::IRegression::run
virtual bool run()
run the computations.
Definition STK_IRegression.h:208
STK::MultidimRegression
The MultidimRegression class allows to regress a multidimensional output variable among a multivariat...
Definition STK_MultidimRegression.h:52
STK
The namespace STK is the main domain space of the Statistical ToolKit project.
main
int main(int argc, char **argv)
Definition tutoAccessors.cpp:3
STK++ QR decomposition:
-----------------------
R matrix:
4.12311 1.21268 5.57832  6.54846
      0 4.41921 2.08981  4.99158
      0       0   4.745  4.22321
      0       0       0 -1.53827
      0       0       0        0
QR matrix:
0  1  2 3
2  3  4 5
2  1  6 7
0  3 -1 2
3 -1  1 1

lapack QR decomposition:
--------------------------
R matrix:
-4.12311 -1.21268 -5.57832 -6.54846
       0 -4.41921 -2.08981 -4.99158
       0        0   -4.745 -4.22321
       0        0        0 -1.53827
       0        0        0        0
QR matrix:
0  1  2 3
2  3  4 5
2  1  6 7
0  3 -1 2
3 -1  1 1
Note
By default the matrix Q is represented as a product of elementary reflectors $ Q = H_1 H_2 . . . H_k$, where $k = \min(m,n)$ each $ H_i $ has the form $ H_i = I - \tau v * v' $. It is possible to get the Q matrix (of size (m,m)) by using the compQ() method. Q will be overwritten.

It is possible to update a QR decomposition when a column is added or removed to the original matrix. In the following example, we remove the second column of matrix and then insert a column with value 1 on the second column.

ExampleOutput
#include <STKpp.h>
using namespace STK;
int main(int argc, char** argv)
{
ArrayXX a(5,4), QR;
a << 0, 1, 2, 3,
2, 3, 4, 5,
2, 1, 6, 7,
0, 3,-1, 2,
3,-1, 1, 1;
stk_cout << "lapack QR decomposition:\n";
lapack::Qr lq(a); lq.run();
// remove column
lq.eraseCol(2);
stk_cout << "R matrix:\n";
stk_cout << lq.R();
QR = lq.Q() * lq.R();
stk_cout << "QR matrix:\n";
stk_cout << QR;
// insert constant column with value 1 in column 2
lq.insertCol(Const::Vector<Real>(5), 2);
stk_cout << "\nR matrix:\n";
stk_cout << lq.R();
QR = lq.Q() * lq.R();
stk_cout << "QR matrix:\n";
stk_cout << QR;
return 0;
}
STK::Const::Vector
Define the constant point.
Definition STK_Const_Arrays.h:601
lapack QR decomposition:
R matrix:
-4.12311 -1.21268 -6.54846
       0 -4.41921 -4.99158
       0        0  4.49464
       0        0        0
       0        0        0
QR matrix:
0  1 3
2  3 5
2  1 7
0  3 2
3 -1 1

R matrix:
-4.12311 -1.21268 -1.69775 -6.54846
       0 -4.41921 -1.11811 -4.99158
       0        0 0.931381  1.39707
       0        0        0   -4.272
       0        0        0        0
QR matrix:
0  1 1 3
2  3 1 5
2  1 1 7
0  3 1 2
3 -1 1 1
Note
the printing mechanism of stk++ does not round small values to zero.

SVD decomposition

The singular-value decomposition of an (m,n) real (or complex) matrix M is a factorization of the form $ \mathbf{U\Sigma V^*}$, where U is an (m,m) real (or complex) unitary matrix, $\mathbf{\Sigma}$ is an (m,n) rectangular diagonal matrix with non-negative real numbers on the diagonal, and $\mathbf{V}$ is an (n,n) real (or complex) unitary matrix.

The diagonal entries $\sigma_i$ of $\mathbf{\Sigma}$ are known as the singular values of M.

ExampleOutput
#include <STKpp.h>
using namespace STK;
int main(int argc, char** argv)
{
ArrayXX a(5,4), usvt;
a << 0, 1, 2, 3,
2, 3, 4, 5,
2, 1, 6, 7,
0, 3,-1, 2,
3,-1, 1, 1;
stk_cout << _T("STK++ Svd decomposition:\n");
stk_cout << _T("------------------------\n");
Svd<ArrayXX> s(a); s.run();
stk_cout << _T("Singular values:\n");
stk_cout << s.D();
stk_cout << _T("\nUDV^T matrix:\n");
stk_cout << s.U()*s.D()*s.V().transpose();
stk_cout << _T("\nlapack Svd decomposition:\n");
stk_cout << _T("---------------------------\n");
lapack::Svd ls(a); ls.run();
stk_cout << _T("Singular values:\n");
stk_cout << ls.D().transpose();
stk_cout << _T("\nUDV^T matrix:\n");
stk_cout << ls.U()*ls.D().diagonalize()*ls.V().transpose();
return 0;
}
STK++ Svd decomposition:
------------------------
Singular values:
12.6179       0       0       0
      0 4.17637       0       0
      0       0 2.51817       0
      0       0       0 1.00223

UDV^T matrix:
 2.5327e-16  1  2 3
          2  3  4 5
          2  1  6 7
5.13478e-16  3 -1 2
          3 -1  1 1

lapack Svd decomposition:
---------------------------
Singular values:
12.6179 4.17637 2.51817 1.00223

UDV^T matrix:
 1.07206e-15  1  2 3
           2  3  4 5
           2  1  6 7
-2.35922e-16  3 -1 2
           3 -1  1 1
Note
Singular values are stored in a vector in lapack method and in a diagonal matrix in STK++ method. It is also possible to compute only U and/or V matrix.

Eigenvalues decomposition

Let A be a square (N,N) matrix with N linearly independent eigenvectors, $ q_i \,\, (i = 1, \dots, N).$ Then A can be factorized as

\[ \mathbf{A}=\mathbf{Q}\mathbf{\Lambda}\mathbf{Q}^{-1} \]

where Q is the square (N,N) matrix whose ith column is the eigenvector $ q_i$ of A and Λ is the diagonal matrix whose diagonal elements are the corresponding eigenvalues, i.e., $ \Lambda_{ii}=\lambda_i $.

If A is a symmetric matrix then Q is an orthogonal matrix.

STK provide native and Lapack interface classes allowing to compute the eigenvalue decomposition of a symmetric square matrix.

ExampleOutput
#include <STKpp.h>
using namespace STK;
int main(int argc, char** argv)
{
CSquareX a(5);
a << 0, 1, 2, 3, 4,
2, 3, 4, 5, 6,
2, 1, 6, 7, 8,
0, 3,-1, 2, 3,
3,-1, 1, 1, 0;
stk_cout << "STK++ Eigen decomposition:\n";
stk_cout << "--------------------------\n";
SymEigen<CSquareX> s(a.upperSymmetrize()); s.run();
stk_cout << "D matrix:\n";
stk_cout << s.eigenValues();
CSquareX QDQt = s.eigenVectors() * s.eigenValues().diagonalize() * s.eigenVectors().transpose();
stk_cout << "QDQ^T matrix:\n";
stk_cout << QDQt;
stk_cout << "\nlapack Eigen decomposition:\n";
stk_cout << "--------------------------\n";
lapack::SymEigen<CSquareX> ls(a); ls.run();
stk_cout << "D matrix:\n";
stk_cout << ls.eigenValues();
QDQt = ls.eigenVectors() * ls.eigenValues().diagonalize() * ls.eigenVectors().transpose();
stk_cout << "QDQ^T matrix:\n";
stk_cout << QDQt;
return 0;
}
STK++ Eigen decomposition:
--------------------------
D matrix:
 20.8996
0.340126
-0.39673
-1.64329
-8.19967
QDQ^T matrix:
-2.33147e-15 1 2 3           4
           1 3 4 5           6
           2 4 6 7           8
           3 5 7 2           3
           4 6 8 3 1.64313e-14

lapack Eigen decomposition:
--------------------------
D matrix:
-8.19967
-1.64329
-0.39673
0.340126
 20.8996
QDQ^T matrix:
-6.88338e-15 1 2 3            4
           1 3 4 5            6
           2 4 6 7            8
           3 5 7 2            3
           4 6 8 3 -7.10543e-15
Note
STK++ eigenvalues computation need a full symmetric matrix as input while lapack version use only upper part of the input data. It is also possible to use the lower part of the matrix in the lapack version.

Solving least square problems

In linear regression, the observations are assumed to be the result of random deviations from an underlying relationship between the dependent variables y and independent variable x.

Given a data set

\[ \{y_{i1},\ldots,y_{id},\, x_{i1}, \ldots, x_{ip}\}_{i=1}^n \]

of n statistical units, a linear regression model assumes that the relationship between the dependent variable y and the regressors x is linear. This relationship is modeled through a disturbance term that adds "noise" to the linear relationship between the dependent variable and regressors. Thus the model takes the form

\[ \mathbf{y}_i = \mathbf{x}^\mathsf{T}_i \boldsymbol\beta + \varepsilon_i, \qquad i = 1, \ldots, n. \]

Often these n equations are stacked together and written in matrix notation as

\[ \mathbf{Y} = X\boldsymbol\beta + \boldsymbol\varepsilon, \, \]

STK++ provide native and Lapack interface classes allowing to solve the least square regression problem.

ExampleOutput
#include <STKpp.h>
using namespace STK;
int main(int argc, char** argv)
{
CArrayXX y(1000,3), x(1000,5), beta(5,3);
Law::Normal l(0,1);
x.rand(l);
beta << 0, 1, 2,
2, 3, 4,
2, 1, 6,
0, 3,-1,
3,-1, 1;
y = x * beta + CArrayXX(1000, 3).rand(l);
stk_cout << "STK++ MultiLeastSquare:\n";
stk_cout << "-----------------------\n";
STK::MultiLeastSquare<CArrayXX, CArrayXX> ols(y,x); ols.run();
stk_cout << "beta matrix:\n";
stk_cout << ols.x();
stk_cout << "\nlapack MultiLeastSquare:\n";
stk_cout << "------------------------\n";
STK::lapack::MultiLeastSquare<CArrayXX, CArrayXX> ls(y,x); ls.run();
stk_cout << "beta matrix:\n";
stk_cout << ls.x();
return 0;
}
STK::ArrayBase::rand
Derived & rand(Law::IUnivLaw< Type > const &law)
set random values to this using a distribution law given by the user.
Definition STK_ArrayBaseApplier.h:122
STK::Law::Normal
Normal distribution law.
Definition STK_Law_Normal.h:72
STK::lapack::MultiLeastSquare
The class MultiLeastSQquare solve the least square problem when the response b is multidimensional.
Definition STK_lapack_MultiLeastSquare.h:102
STK++ MultiLeastSquare:
-----------------------
beta matrix:
  0.0247875   1.02139  1.99183
    1.96609   3.01792   4.0079
    2.04481   1.03293  6.03462
-0.00848336   2.99758 -1.05867
    3.02074 -0.996975  1.00954

lapack MultiLeastSquare:
------------------------
beta matrix:
  0.0247875   1.02139  1.99183
    1.96609   3.01792   4.0079
    2.04481   1.03293  6.03462
-0.00848336   2.99758 -1.05867
    3.02074 -0.996975  1.00954
Note
It is also possible to solve weighted least-square regression problems.

Inverting matrices

Matrices can be inverted using either the templated functor STK::InvertMatrix or the templated function STK::invert. The first example below deals with general square and/or symmetric matrices.

ExampleOutput
#include <STKpp.h>
using namespace STK;
int main(int argc, char** argv)
{
CArray<Real, 4, 4> a(4,4);
a << 0, 1, 2, 3,
2, 3, 4, 5,
2, 1, 6, 7,
0, 3,-1, 2;
stk_cout << "Inverse general matrix:\n";
stk_cout << "-----------------------\n";
stk_cout << a*invert(a);
stk_cout << "\nInverse upper-symmetric matrix:\n";
stk_cout << "-------------------------------\n";
stk_cout << a.upperSymmetrize();
stk_cout << a.upperSymmetrize()*invert(a.upperSymmetrize());
stk_cout << "\nInverse lower-symmetric matrix:\n";
stk_cout << "-------------------------------\n";
stk_cout << a.lowerSymmetrize();
stk_cout << a.lowerSymmetrize()*invert(a.lowerSymmetrize());
return 0;
}
Inverse general matrix:
-----------------------
           1 0 0 0
           0 1 0 0
-8.88178e-16 0 1 0
           0 0 0 1

Inverse upper-symmetric matrix:
-------------------------------
0 1 2 3
1 3 4 5
2 4 6 7
3 5 7 2
           1 0            0            0
           0 1 -2.22045e-16 -2.22045e-16
 6.66134e-16 0            1 -2.22045e-16
-2.10942e-15 0  1.11022e-16            1

Inverse lower-symmetric matrix:
-------------------------------
0 2  2  0
2 3  1  3
2 1  6 -1
0 3 -1  2
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1

The second example deals with lower and upper triangular arrays.

ExampleOutput
#include <STKpp.h>
using namespace STK;
int main(int argc, char** argv)
{
Array2DLowerTriangular<Real> a(5,5);
a << 1,
1, 2,
3, 4, 3,
4, 5, 6, 6,
7, 8, 2, 3, 2;
stk_cout << "Inverse lower-triangular matrix:\n";
stk_cout << "--------------------------------\n";
stk_cout << a*invert(a);
stk_cout << "\nInverse upper-triangular matrix:\n";
stk_cout << "--------------------------------\n";
stk_cout << a.transpose()*invert(a.transpose());
return 0;
}
Inverse lower-triangular matrix:
--------------------------------
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0
0 0 0 0 1

Inverse upper-triangular matrix:
--------------------------------
1 0 0 -1.11022e-16 0
0 1 0            0 0
0 0 1            0 0
0 0 0            1 0
0 0 0            0 1
Note
If the matrix is not invertible, the result provided will be a generalized inverse.