% [f,fX,E] = rrrtrain(X[,Y,L,l,m]) Train RRR linear function f(x) = W1.W2.x+w
%
% Reduced-rank regression (RRR) is a least-squares linear regression with a
% low-rank parameter matrix:
%   min_{W2,W1,w}{ sum_n{ |yn - f(xn)|² } }
% where f(x) = W1.W2.x+w with W1 of D1xL and W2 of LxD2, with L <= D1.
% Hence the regression matrix C = W1.W2 is low-rank. Storing and operating
% with (W1,W2) is more efficient than with C. If X=Y RRR becomes PCA.
%
% Solution for RRR: let XX, YY and XY be the covariance matrices of (x,x),
% (y,y) and (x,y), and the RRR matrix be RRR = XY'.inv(XX).XY. Then the
% full-rank (regular) regression matrix is W = XY'.inv(XX) and the RRR
% regression matrices are W1 = U(:,1:L)), W2 = W1'.W where (U,s) = svd(RRR)
% are the eigenvectors and eigenvalues of RRR (in decreasing order). So, given
% U and W we can construct the RRR solution for any rank L.
% If using a regularisation parameter l>0, we use XX+l.I (which is always
% full-rank) instead of XX.
%
% In:
%   X: NxD2 matrix, N D2-dim data points rowwise, the training input vectors.
%   Y: NxD1 matrix, N D1-dim data points rowwise, the training output vectors.
%      When X=Y, use Y=[] (faster).
%   L: desired rank, should be <= D1.
%   l: (nonnegative scalar) regularisation parameter. Default: 0.
%   m: if 1, return in f also the full-rank regression matrix W and bias ww,
%      the eigenvectors U and eigenvalues s of the RRR matrix, and the fit
%      error e for L=1:D1 (useful to select the value of L). Default: 0.
% Out:
%   f: (struct) the RRR linear function, with fields:
%      type='rrr', W1 (D1xL), W2 (LxD2), w (D1x1), regularisation parameter l.
%      Optionally: W (D1xD2), ww (D1x1), U (D1xD1), s (D1x1), e (D1x1).
%   fX: NxD1 matrix, f(X).
%   E: fit error.
%
% Any non-mandatory argument can be given the value [] to force it to take
% its default value.

% Copyright (c) 2016 by Miguel A. Carreira-Perpinan

function [f,fX,E] = rrrtrain(X,Y,L,l,m)

% ---------- Argument defaults ----------
if ~exist('l','var') l = []; f.l = 0; else f.l = l; end;
if ~exist('m','var') || isempty(m) m = 0; end;
% ---------- End of "argument defaults" ----------

[N,D2] = size(X); f.type = 'rrr';

X1 = sum(X,1)'; XX = X'*X - X1*(X1'/N);
if ~isempty(l) XX = XX + spdiags(repmat(l,D2,1),0,D2,D2); end
if ~isempty(Y)
  Y1 = sum(Y,1)'; YY = Y'*Y - Y1*(Y1'/N); XY = X'*Y - X1*(Y1'/N);
  W = XY' / XX; RRR = W*XY; RRR = (RRR+RRR')/2; [U,s] = svd(RRR);
  f.W1 = U(:,1:L); f.W2 = f.W1'*W;
  if m
    f.W = W; f.U = U; f.s = diag(s)/N;
    f.ww = (Y1-f.W*X1)/N;		% = mean(Y-X*f.W',1)'
    % compute fit error for every L=1:D1
    Y0 = bsxfun(@minus,Y,Y1'/N); f.e = sum(sum(Y0.^2)); Y0 = Y0*U;
    WX0 = bsxfun(@minus,X*f.W',(X1'*f.W')/N)*U;
    f.e = f.e + cumsum(sum(WX0.^2,1) - 2*sum(Y0.*WX0,1))';
  end;
else		% if Y = X then RRR = PCA; the following is more efficient
  [U,s] = svd((XX+XX')/2); f.W1 = U(:,1:L); f.W2 = f.W1';
  if m
    f.W = eye(D2,D2); f.ww = zeros(D2,1); f.U = U; f.s = diag(s)/N;
    f.e = N*(sum(f.s) - cumsum(f.s));
  end;
  Y1 = X1;
end
f.w = (Y1-(f.W1*(f.W2*X1)))/N;		% = mean(Y-((X*f.W2')*f.W1'),1)'

if nargout>1 fX = rrr(X,f); end
if nargout>2 E = sum(sum((Y-fX).^2)); end

