function [Glambda,lambda,X] = invregfm(U,s,V,Y,lambda);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [Glambda,lambda,X] = invreg(U,s,V,Y,lambda);
% finds optimal regularization parameter lambda 
%
%      lambda =   arg min G(l)= || P(l) Y ||^2 ./ Trace(n P(l))^2 
%
% for using in Tikhonov regularization
%
%         min { || A x - Y ||^2 + lambda^2 ||L x ||^2 } (function tikh), 
%
% where        
%
%      P(l) = I-A inv(A' A+l^2 I) A'
%
% and
%
%      A = U diag(s) V' by compact SVD (csvd). 
%
% The seudoinverse of A' A, A' A+lambda^2 I, is implicitly used.
%
% Input arguments:
%
%   U      - left orthognal matrix from SVD of A
%   s      - column vector of singular values of A
%   V      - right orthognal matrix from SVD of A (if only the first r columns of V are passed
%          - then only the first r components of x are calculated)
%   Y      - data matrix
%   lambda - regularization parameter when it is input
%
% Output arguments:
%
%   Glambda - G(lambda)
%   lambda  - optimal regularization parameter
%   X       - solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%Comprobar si los autovalores estan ordenados, en otro caso da error
%%% WARNING: Usa fmin

[sa,i]=sort(s);
if ~all(i == [length(i):-1:1]')
  disp('ordenando');
  s=sa;
  s=flipud(s(:));
  i=flipud(i(:));
  U=U(:,i);
  V=V(:,i);
end

umbral=1e-8;
[n,m]=size(Y);
p=length(s);
UtY=U'*Y;
UtYc=Y-U*UtY;
onesm=ones(m,1);
t=abs(UtY).^2*onesm;
tc=abs(UtYc).^2*onesm;
s2=s.^2;
tolerance=s2(1)*umbral;
indexes=find(s2>=tolerance);
matrixrank=length(indexes);

rankrange=1:matrixrank;
UtY=UtY(rankrange,:);
s2=s2(rankrange);
sumt=sum(t(matrixrank+1:p)) + sum(tc);
t=t(rankrange);
if nargin<5
  npoints=100;
  minregpar=s(1)*eps;
  maxregpar=minregpar.*((1.2^npoints)/eps);
  opt=foptions;
  opt(1) = 0;
  opt(2) = 1.0e-10;
  opt(3) = 1.0e-10;

%%tgcvgen

  lambda = fmin('gcvgen',log10(minregpar),log10(maxregpar),opt,n,s2,t,sumt);
  lambda = 10.^lambda;
end;
  
Glambda=gcvgen(log10(lambda),n,s2,t,sumt);
if nargout>2
  X=V(:,rankrange)*(((s(rankrange)./(s2+lambda.^2))*onesm').*UtY);
end;
