function [B, S, h, peff] = mridge(Cxx, Cyy, Cxy, dof, options);%MRIDGE  Multiple ridge regression with generalized cross-validation.%%   [B, S, h, peff] = MRIDGE(Cxx, Cyy, Cxy, dof, OPTIONS) returns a%   regularized estimate B = Mxx_h Cxy of the coefficient matrix in%   the multivariate multiple regression model Y = X*B + noise(S). The%   matrix Mxx_h is the regularized inverse of the covariance matrix%   Cxx,%%             Mxx_h = inv(Cxx + h^2 * I).%%   The matrix Cxx is an estimate of the covariance matrix of the%   independent variables X, Cyy is an estimate of the covariance%   matrix of the dependent variables Y, and Cxy is an estimate of the%   cross-covariance matrix of the independent variables X and the%   dependent variables Y. The scalar dof is the number of degrees of%   freedom that were available for the estimation of the covariance%   matrices.%%   The input structure OPTIONS contains optional parameters for the%   algorithm:%%     Field name         Parameter                                   Default%%     OPTIONS.regpar     Regularization parameter h. If regpar       not set%                        is set, the scalar OPTIONS.regpar is%                        taken as the regularization parameter h. %                        If OPTIONS.regpar is not set (default), %                        the regularization parameter h is selected %                        as the minimizer of the generalized %                        cross-validation (GCV) function. The output%                        variable h then contains the selected %                        regularization parameter.%  %     OPTIONS.relvar_res Minimum relative variance of residuals.       5e-2%                        From the parameter OPTIONS.relvar_res, a%                        lower bound for the regularization parameter%                        is constructed, in order to prevent GCV from%                        erroneously choosing too small a %                        regularization parameter (see GCVRIDGE).%%   The OPTIONS structure is also passed to GCVRIDGE.%%   MRIDGE returns the ridge estimate B of the matrix of regression%   coefficients. Also returned are an estimate S of the residual%   covariance matrix, the regularization parameter h, and the scalar%   peff, an estimate of the effective number of adjustable%   parameters in B.  %  %   MRIDGE computes the estimates of the coefficient matrix and of the%   residual covariance matrix from the covariance matrices Cxx, Cyy,%   and Cxy by solving the regularized normal equations. The normal%   equations are solved via an eigendecomposition of the covariance%   matrix Cxx. However, if the data matrices X and Y are directly%   available, a method based on a direct factorization of the data%   matrices will usually be more efficient and more accurate.%%   See also: IRIDGE, GCVRIDGE.  error(nargchk(4, 5, nargin))     % check number of input arguments     px           = size(Cxx, 1);  py           = size(Cyy, 1);  if size(Cxx, 2) ~= px | size(Cyy, 2) ~= py | any(size(Cxy) ~= [px, py])     error('Incompatible sizes of covariance matrices.')  end    % ==============           process options        ========================  if nargin == 4 | isempty(options)    options    = [];    fopts      = [];  else    fopts      = fieldnames(options);  end      if strmatch('regpar', fopts)    h          = options.regpar;     h_given    = 1;  else    h_given    = 0;  end    if strmatch('relvar_res', fopts)    relvar_res = options.relvar_res;   else    relvar_res = 5e-2;  end  % =================           end options        =========================    if nargout > 1    S_out      = 1;  else    S_out      = 0;  end      % eigendecomposition of Cxx  rmax         = min(dof, px);     % maximum possible rank of Cxx  [V, d, r]    = peigs(Cxx, rmax);    % Fourier coefficients. (The following expression for the Fourier  % coefficients is only correct if Cxx = X'*X and Cxy = X'*Y for  % some, possibly scaled and augmented, data matrices X and Y; for  % general Cxx and Cxy, all eigenvectors V of Cxx must be included,  % not just those belonging to nonzero eigenvalues.)  F            = repmat(ones(r, 1)./sqrt(d), 1, px) .* V' * Cxy;  % Part of residual covariance matrix that does not depend on the  % regularization parameter h:  if (dof > r)     S0         = Cyy - F'*F;  else    S0         = sparse(py, py);  end    if ~h_given    % approximate minimum squared residual    trSmin     = relvar_res * trace(Cyy);        % find regularization parameter that minimizes the GCV object function    h          = gcvridge(F, d, trace(S0), dof, r, trSmin, options);  end    % get matrix of regression coefficients  B            = V * (repmat(sqrt(d)./(d + h^2), 1, py) .* F);    if S_out    % get estimate of covariance matrix of residuals    S          = S0 + F' * (repmat(h^4./(d + h^2).^2, 1, py) .* F);  end    if nargout == 4    % effective number of adjusted parameters: peff = trace(Mxx_h Cxx)    peff       = sum(d ./ (d + h^2));  end