Factorizations: LU, QR, Cholesky, singular values, eigenvalues, Hessenberg

Functions

array lu (const array &in)
 LU factorization (packed). More...
 
void lu (array &lower, array &upper, const array &in)
 LU factorization. More...
 
void lu (array &lower, array &upper, array &pivot, const array &in)
 LU factorization (with pivoting). More...
 
array qr (const array &in)
 QR factorization (packed). More...
 
void qr (array &q, array &r, const array &in)
 QR factorization. More...
 
void qr (array &q, array &r, array &tau, const array &in)
 QR factorization with tau. More...
 
array cholesky (const array &X, bool is_upper=true)
 Cholesky decomposition ("Y^T * Y == X"). More...
 
array cholesky (unsigned &info, const array &X, bool is_upper=true)
 Cholesky decomposition ("Y^T * Y == X"). More...
 
array hessenberg (const array &in)
 Hessenberg matrix form. More...
 
void hessenberg (array &h, array &q, const array &in)
 Hessenberg matrix h with unitary permutation matrix q. More...
 
array eigen (const array &in, bool is_diag=false)
 Compute eigenvalues. More...
 
void eigen (array &values, array &vectors, const array &in)
 Eigenvalues and eigenvectors. More...
 
array svd (const array &in, bool is_diag=false)
 Singular values. More...
 
void svd (array &s, array &u, array &v, const array &in)
 Singular values with unitary bases: in = u * s * v. More...
 

Detailed Description

Function Documentation

array af::lu ( const array &  in)

LU factorization (packed).

Double-precision or complex values require ArrayFire Pro.

Parameters
[in]in
Returns
lower and upper matrix combined (getrf)
Examples:
examples/getting_started/lin_algebra.cpp.
void af::lu ( array &  lower,
array &  upper,
const array &  in 
)

LU factorization.

Double-precision or complex values require ArrayFire Pro.

array l, u;
array in = randu(5, 5);
lu(l, u, in); // Unpacked output
print(l);
print(u);
Parameters
[out]lowertriangular matrix
[out]uppertriangular matrix
[in]in
void af::lu ( array &  lower,
array &  upper,
array &  pivot,
const array &  in 
)

LU factorization (with pivoting).

This function decomposes a matrix A as a product of a lower-triangular matrix l and an upper-triangular matrix u. Note that in general, the product of l and u gives us a matrix which is a row-permutation of the original matrix A, obtained by pre-multiplying A with a permutation matrix P.

The following relationship is obtained as a result:

lu = PA

The permutation matrix P is simply a row-permutation of the identity matrix. The location of non-zero entries (1's) in the P matrix can be obtained via a permutation vector p returned by this function, as explained in the code example below.

Double-precision or complex values require ArrayFire Pro.

#include<stdio.h>
#include <arrayfire.h>
using namespace af;
int main()
{
// create a lower triangular matrix L from a host buffer
float h_buff[] = {1, 1, 2, 0, 4, 5, 0, 0, 8};
array L(3, 3, h_buff); // L = 1 0 0
print(L); // 1 4 0
// 2 5 8
// create an upper triangular matrix U from a host buffer
float g_buff[] = {7, 0, 0, 3, 4, 0, 6, 7, 8};
array U(3, 3, g_buff); // U = 7 3 6
print(U); // 0 4 7
// 0 0 8
// matrix multiplication
array C = matmul(L, U); // C = LU = 7 3 6
print(C); // 7 19 34
// 14 26 111
// reverse engineer l and u from C
array l, u, p;
lu(l, u, p, C);
print(l); // lower triangular matrix with diagonals = 1
print(u); // upper triangular matrix
print(p); // row location of 1's for each column in the permutation matrix
// l = 1.0000 0.0000 0.0000
// 0.5000 1.0000 0.0000
// 0.5000 -0.6000 1.0000
// u = 14.0000 26.0000 111.0000
// 0.0000 -10.0000 -49.5000
// 0.0000 0.0000 -51.2000
// p = 2 0 1
// create permutation matrix based on the indices in p
//
// row row row
// 2 0 1
//
// | | |
// | | |
//
float permutation_matrix[] = { 0, 1, 0,
0, 0, 1,
1, 0, 0 };
array P(3, 3, permutation_matrix);
print(P);
// check if lu = PC
array lu = matmul(l, u);
print(lu);
// lu = 14.0000 26.0000 111.0000
// 7.0000 3.0000 6.0000
// 7.0000 19.0000 34.0000
array PC = matmul(P, C);
print(PC);
// PC = 14.0000 26.0000 111.0000
// 7.0000 3.0000 6.0000
// 7.0000 19.0000 34.0000
return 0;
}
Parameters
[out]lowertriangular matrix
[out]uppertriangular matrix
[out]pivotindices
[in]in
array af::qr ( const array &  in)

QR factorization (packed).

Double-precision or complex values require ArrayFire Pro.

array in = randu(5, 5);
array out = qr(in);
print(out);
Parameters
[in]in
Returns
Q,R factors (geqrf)
void af::qr ( array &  q,
array &  r,
const array &  in 
)

QR factorization.

Double-precision or complex values require ArrayFire Pro.

array q, r;
array in = randu(5, 5);
qr(q, r, in); // Unpacked output
print(q);
print(r);
Parameters
[out]qunitary rotation matrix
[out]rupper triangular
[in]in
void af::qr ( array &  q,
array &  r,
array &  tau,
const array &  in 
)

QR factorization with tau.

Double-precision or complex values require ArrayFire Pro.

Parameters
[out]qorthogonal matrix
[out]rupper triangular matrix
[out]tauAdditional information about q
[in]in
array af::cholesky ( const array &  X,
bool  is_upper = true 
)

Cholesky decomposition ("Y^T * Y == X").

Double-precision or complex values require ArrayFire Pro.

float h_x[]={ 123.1029, -15.4097, -15.4097, 62.9234 };
array x = array(2,2,h_x);
array y = cholesky(x,true);
print(x);
print(y);
// x =
// 123.1029 -15.4097
// -15.4097 62.9234
// y =
// 11.0952 -1.3889
// 0.0000 7.8099
Parameters
[in]X
[in]is_uppertrue for upper triangular solution, false for lower triangular solution
Returns
upper or lower triangular factor
Note
Exception is throw in case input is not positive definite.
array af::cholesky ( unsigned &  info,
const array &  X,
bool  is_upper = true 
)

Cholesky decomposition ("Y^T * Y == X").

Double-precision or complex values require ArrayFire Pro.

Parameters
[out]infodetails on result of decomposition. If successful info is 0 indicating input is positive definite. Else a positive value representing the row at which decomposition failed.
[in]X
[in]is_uppertrue for upper triangular solution, false for lower triangular solution
Returns
upper or lower triangular factor
Note
No exception is thrown for this entry point.
array af::hessenberg ( const array &  in)

Hessenberg matrix form.

Double-precision or complex values require ArrayFire Pro.

Parameters
[in]in
Returns
Hessenberg form
void af::hessenberg ( array &  h,
array &  q,
const array &  in 
)

Hessenberg matrix h with unitary permutation matrix q.

such that in = q * h * q^T. Double-precision or complex values require ArrayFire Pro.

Parameters
[out]hHessenberg matrix
[out]qunitary matrix
[in]in
array af::eigen ( const array &  in,
bool  is_diag = false 
)

Compute eigenvalues.

Requires ArrayFire Pro.

Parameters
[in]in
[in]is_diagtrue then produce diagonal matrix of eigenvalues, false then vector of eigenvalues
Returns
eigenvalues
Examples:
examples/getting_started/lin_algebra.cpp.
void af::eigen ( array &  values,
array &  vectors,
const array &  in 
)

Eigenvalues and eigenvectors.

Requires ArrayFire Pro.

This function returns eigenvalues and eigenvectors for the input matrix

#include<stdio.h>
#include <arrayfire.h>
using namespace af;
int main(void)
{
try {
// create a coefficient matrix A
float h_buffer[] = {1, -4, 6, -4, 2, -2, 6, -2, 1};
array A(3, 3, h_buffer);
print(A);
// A = 1 -4 6
// -4 2 -2
// 6 -2 1
// calculate eigenvalues and eigenvectors
array values, vectors;
eigen(values, vectors, A);
// In this case, the three (scalar) Eigenvalues
// are located on the diagonal
print(values); // values = -5.3401 0.0000 0.0000
// 0.0000 -0.1188 0.0000
// 0.0000 0.0000 9.4588
// The Eigenvectors are returned as column vectors
// Each vector corresponds to the Eigenvalue in that
// particular column
print(vectors); // vectors = 0.7420 0.1805 -0.6457
// 0.2330 0.8336 0.5008
// -0.6286 0.5221 -0.5764
} catch(af::exception &ae) {
std::cout << ae.what() << std::endl;
}
return 0;
}

In the above code example, we obtain three distinct eigenvalues and three distinct eigenvectors, as illustrated in the figure below. As a check, if we pre-multiply each eigenvector with the input matrix, the effect is identical to scaling the eigenvector by a non-zero scalar. This scalar is the associated eigenvalue.

eigen.png
Parameters
[out]values
[out]vectors
[in]in
array af::svd ( const array &  in,
bool  is_diag = false 
)

Singular values.

Double-precision or complex values require ArrayFire Pro.

array in = randu(5, 5);
array s = svd(in);
print(s);
// s =
// 2.6075
// 0.9519
// 0.4444
// 0.1295
// 0.1070
Parameters
[in]in
[in]is_diagtrue then produce diagonal matrix of eigenvalues, false then vector of eigenvalues
Returns
singular values
void af::svd ( array &  s,
array &  u,
array &  v,
const array &  in 
)

Singular values with unitary bases: in = u * s * v.

Double-precision or complex values require ArrayFire Pro.

array s, u, v;
array in = randu(5, 5);
svd(s, u, v, in); // Unpacked output
print(s); // 5x5 diag matrix
print(u);
print(v);
Parameters
[out]ssingular values
[out]uleft unitary matrix
[out]vright unitary matrix
[in]inThe input matrix.