Interpolative Decomposition based on Strong RRQR

A matlab implementation of the interpolative decomposition based on strong RRQR

The interpolative decomposition (ID) (see reference in the bottom) is a useful low-rank form in many structured matrix representations. This matlab code has been used a lot in my research for an efficient ID approximation.

Given an m*n matrix A, the code gives an rank-k ID approximation of A in the form of

A \approx U * A(J, :)

where J is a row index subset of size k, and U is a m*k matrix with entries bounded by the pre-specified constant. A(J,:) and U are usually called the skeleton and projection matrix, respectively.

The code contains the ID approximation with (1) a fixed rank, i.e., the dimension of R11; (2) an error threshold, i.e., each row error vector of the above approximation has norm less than the threshold.

ALGORITHM DESCRIPTION:

Denote the partial strong RRQR of A' as

A' * P = [Q1, Q2] * [R11, R12; 0, R22]

\approx Q1 * [R11, R12]

where P is a permutation matrix, [Q1, Q2] is an orthgonal matrix, R11 is an upper-triangular matrix, and R22 is not necessarily an upper-triangular matrix. In addition, by strong RRQR, entries of the matrix inv(R11)*R12 all have absolute value bounded by a pre-specified constant that is not less than 1.

The ID approximation of A can be obtained as

A' * P \approx Q1 * [R11, R12]

= [Q1*R11, Q1*R12]

= [A1', A1'*inv(R11)*R12] (A1' denotes the corresponding first few columns of A'P)

= A1' * [eye, inv(R11)*R12]

= A1' * [eye, E'] (denote inv(R11)*R12 as E')

Then, by simple transformation

P' * A \approx [eye; E] * A1

A \approx P * [eye; E] *A1

Denote U = P*[eye;E] and A1 = A(J, :), the ID is written as

A \approx U * A(J, :);

Reference:

Cheng, Hongwei, Zydrunas Gimbutas, Per-Gunnar Martinsson, and Vladimir Rokhlin. "On the compression of low rank matrices." SIAM Journal on Scientific Computing 26, no. 4 (2005): 1389-1404.

Note:

QRpartial.c is a simplified partial pivoting QR code which do not return the orthogonal matrix Q and only returns the permutation vector and the half-finished triangular matrix [R11, R12; 0, R22]. (R22 is not necessary to be an upper-triangular matrix).

QRpartial.c can be faster than full qr.m in Matlab if rank of ID is small. However, it still cannot compete with the full matlab qr.m in some cases. There is still chance to further accelerate QRpartial.c simply through a better coding. Greatly appreciated if any one can provide an efficient partial pivoting QR code for both a fixed rank or a fixed threshold.