function [C, G, R] = SkeletonApproximation(A, r)
%
% This function provides a numerical implementation of the Skeleton
% Approximation Algorithm described in [1].
%
% Input:
% - A: m x n matrix.
% - r: number of columns and rows to be selected from A.
%
% Ouput:
% - C: m x r matrix containing r columns from A.
% - R: r x n matrix containing r rows from A.
% - G: r x r matrix.
%
% Outcome:
% The rank-r matrix T = CGR is ''close'' to the input matrix A, i.e. the
% approximation error ||A - T||_2 is ''close'' to ||A - Ar||_2, where Ar
% is the ''best'' rank-r matrix computed by the Singular Value
% Decomposition of A.
%
% Notes:
% The authors of [1] give Theorem 3.1 and describe a constructive proof
% for their Theorem. The implementation that we describe here is a
% special case of Theorem 3.1; in particular we set F = A - Ar, and the
% matrices \hat{U} and \hat{V} of the Theorem are constructed by running
% the pivoted QR decomposition method of [2]. This column/row selection
% step can be replaced by any Rank Revealing QR algorithm with similar
% properties (see, for example, Table 2 in [3]).
%
% References:
%
% [1] S. A. Goreinov, E. E. Tyrtyshnikov, and N. L. Zamarashkin.
% A theory of pseudoskeleton approximations.
% Linear Algebra and its applications, Volume 261, Issues 1-3, August
% 1997, Pages 1-21
%
% [2] Gene Golub.
% Numerical methods for solving linear least squares problems.
% Numerische Mathematik, 7:206–216, 1965.
%
% [3] C. Boutsidis, M. W. Mahoney, and P. Drineas.
% An improved approximation algorithm for the column subset selection
% problem.
% Proceedings of the Twentieth Annual ACM-SIAM Symposium on Discrete
% Algorithms, SODA 2009, New York, NY, USA, January 4-6, 2009.
%
% -------------------------------------------------------------------------
% Author: Christos Boutsidis, August 2009.
% E-mail: christos.boutsidis@gmail.com
% -------------------------------------------------------------------------
% The first step of the algorithm (see eqn. 3.2 in [1]) requires the
% computation of the Singular Value Decomposition of the matrix A - F.
% Since F = A - Ar, it is A - F = Ar, thus we only need to compute the
% top-r singular values and vectors of the matrix A.
[U, S, V] = svds(A,r);
% Construct the matrix C containing r columns from A.
[qV, rV, pV] = qr(V', 0);
C = A(:, pV(1:r));
% Construct the matrix R containing r rows from A.
[qU, rU, pU] = qr(U', 0);
R = A(pU(1:r),:);
% Construct the matrix G. We set G as the solution of the following problem
% min_G || A - CGR ||_F.
% This choice of G is in general different from the one in [1]. We chose
% this G since we find its construction simpler than the one in [1] and we
% believe that - in practice - is as ''good'' as the one in [1].
G = pinv(C, .05) * A * pinv(R, .05);