lowrank_only_pareto.m

function out_lowrank = lowrank_only_pareto(A, nnz_targets)
%LOWRANK_ONLY_PARETO  Low-rank-only baseline: best rank-k approximation of A.
%
%   out_lowrank = LOWRANK_ONLY_PARETO(A) approximates A by its optimal rank-k
%   truncated SVD at a sequence of storage budgets (relative NNZ levels). Each
%   level yields one (relative_NNZ, relative_error) point, so the result is the
%   "low-rank-only" Pareto curve: the accuracy obtainable from a pure low-rank
%   approximation alone, with no sparse part. Comparing it against R_TRUST /
%   R3_TRUST shows how much the sparse component improves accuracy at the same
%   storage cost.
%
%   A rank-r factor pair (H_{m x r}, W_{r x n}) stores r*(m+n) entries. Given a
%   total budget of B entries, the affordable rank is r = floor(B / (m+n)). So at
%   the same relative NNZ as a sparse matrix with B nonzeros, the best low-rank
%   approximation uses rank r = floor(B / (m+n)).
%
%   out_lowrank = LOWRANK_ONLY_PARETO(A, nnz_targets) specifies the budgets as
%   fractions of numel(A) in (0,1]. Default: [0.015 0.020 0.025 0.030].
%
%   The returned struct array mirrors the R_TRUST / R3_TRUST output fields
%   (relative_NNZ, relative_error, time, and additionally rank), so it plugs
%   straight into PLOT_PARETO:
%
%       out_lowrank = lowrank_only_pareto(A);
%       plot_pareto({out_lowrank, out_RTrust}, ...
%                   {'Low-Rank-Only','Sparse+Low-Rank'}, 'YScale','log');
%
%   The approximation error is computed via the Frobenius-norm identity
%   (Lemma 4.1):  ||A - L_k||_F = sqrt( ||A||_F^2 - sum_{i=1}^{k} sigma_i^2 ),
%   which avoids forming the m-by-n product and is exact.

if nargin < 2 || isempty(nnz_targets)
    nnz_targets = [0.005 0.01 0.015 0.020 0.025 0.030];          % 1.5% ... 3.0% of entries
end

[m, n]     = size(A);
normA      = norm(A, "fro");
normA_sq   = normA^2;
nnz_targets = sort(nnz_targets(:)');                   % ascending budget

out_lowrank = struct([]);
prev_rank   = 0;                                       % track last computed rank
j           = 0;                                       % output index (skips duplicates)
for i = 1 : numel(nnz_targets)
    budget  = round(nnz_targets(i) * m * n);           % total entries available
    r       = floor(budget / (m + n));                  % affordable rank
    r       = max(r, 1);                                % at least rank 1
    r       = min(r, min(m, n));                        % at most full rank

    if r == prev_rank                                  % same rank as previous budget
        continue;                                      % skip (result would be identical)
    end
    prev_rank = r;
    j = j + 1;

    fprintf('Low-rank-only %d/%d: budget %.2f%% -> rank %d\n', ...
        i, numel(nnz_targets), 100*nnz_targets(i), r);

    timer = tic;
    [~, S_svd, ~] = svdsecon(A, r);                    % thin SVD (only top-r)
    elapsed = toc(timer);

    sigma_sq_sum = sum(diag(S_svd).^2);                % sum of squared sing. vals

    % Lemma 4.1: ||A - L_k||_F^2 = ||A||_F^2 - sum_{i=1}^k sigma_i^2
    residual_sq  = max(normA_sq - sigma_sq_sum, 0);    % guard against round-off

    out_lowrank(j).rank           = r;
    out_lowrank(j).relative_NNZ   = r * (m + n) / (m * n);
    out_lowrank(j).relative_error = sqrt(residual_sq) / normA;
    %out_lowrank(j).time           = elapsed;
end
end

% ========================================================================
% Fast thin SVD for a small number of singular values
% ========================================================================
function [U, S, V] = svdsecon(X, k)
%SVDSECON  Equivalent to svds(X,k) but faster when k << min(size(X)).
%   Requires k < min(m,n). Works on the smaller of X*X' / X'*X via eigs.
%   Vipin Vijayan (2014).
[m, n] = size(X);
assert(k <= m && k <= n, 'k must be smaller than both dimensions of X');

if m <= n
    C = X * X';
    [U, D] = eigs(C, k);
    clear C;
    if nargout > 2
        V = X' * U;
        s = sqrt(abs(diag(D)));
        V = bsxfun(@(x, c) x ./ c, V, s');
        S = diag(s);
    end
else
    C = X' * X;
    [V, D] = eigs(C, k);
    clear C;
    U = X * V;
    s = sqrt(abs(diag(D)));
    U = bsxfun(@(x, c) x ./ c, U, s');
    S = diag(s);
end
end