Skip to content

MATLAB bindings for CQRRPT and CQRRT #8

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mmelnich wants to merge 1 commit into
base: main
Choose a base branch
from
Open

MATLAB bindings for CQRRPT and CQRRT #8

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion CMakeLists.txt
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -38,5 +38,5 @@ else()
"Set -DMatlab_ROOT_DIR=/path/to/MATLAB to enable. "
"The +randlapack package can still be used by adding "
"${CMAKE_CURRENT_SOURCE_DIR}/matlab to the MATLAB path, "
"but calls into bqrrp_mex will fail until the MEX is built.")
"but calls into the MEX files will fail until they are built.")
endif()
31 changes: 27 additions & 4 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,13 @@ MATLAB, and eventually Python, bindings for [RandLAPACK](https://github.com/Ball

## Scope (v0)

Proof-of-concept binding for a single driver, `BQRRP` (randomized blocked QR with column pivoting), in single and double precision, real matrices, column-major. More drivers as their C++ APIs stabilize. Python next.
MATLAB bindings for three QR drivers, in single and double precision, real matrices, column-major:

* `BQRRP`: randomized blocked QR with column pivoting, any shape.
* `CQRRPT`: randomized rank-revealing QR with column pivoting, for tall matrices.
* `CQRRT`: randomized Cholesky QR without pivoting, for tall full-rank matrices.

More drivers as their C++ APIs stabilize. Python next.

## Requirements

Expand All @@ -24,7 +30,7 @@ cmake --build build -j

Pass `-DCMAKE_BUILD_TYPE=Release`: without a build type the MEX compiles without optimization (Ninja with MSVC even defaults to Debug).

If MATLAB is not found, MEX targets are skipped with a clear warning; the `.m` files still ship and can be used once the MEX is built. On success, `matlab/+randlapack/private/bqrrp_mex.<ext>` is in place.
If MATLAB is not found, MEX targets are skipped with a clear warning; the `.m` files still ship and can be used once the MEX is built. On success, `bqrrp_mex.<ext>`, `cqrrpt_mex.<ext>` and `cqrrt_mex.<ext>` are in `matlab/+randlapack/private/`.

### Windows

Expand Down Expand Up @@ -56,12 +62,29 @@ err = norm(A(:, J) - Q*R, 'fro') / norm(A, 'fro');
Two output modes, chosen by a trailing string (qr-style):

```matlab
[Q, R, J] = randlapack.bqrrp(A, 'explicit') % default; matches qr(A, 'vector')
[Q, R, J] = randlapack.bqrrp(A, 'explicit') % default; matches qr(A, 'econ', 'vector')
[A_out, tau, J] = randlapack.bqrrp(A, 'implicit') % GEQP3-format; skips Q materialization
```

Run `help randlapack.bqrrp` for the full signature, including the optional `b_sz`, `d_factor`, and RNG `state` arguments.

For **tall** matrices (`m >= d_factor*n`, with `d_factor = 1.25` by default), two sketch-and-precondition drivers:

```matlab
A = randn(20000, 200);

% CQRRPT: rank-revealing, pivoted. The detected rank is size(R, 1).
[Q, R, J] = randlapack.cqrrpt(A);
err = norm(A(:, J) - Q*R, 'fro') / norm(A, 'fro');

% CQRRT: unpivoted, for full-rank A. Same layout as qr(A, 0).
[Q, R] = randlapack.cqrrt(A);
err = norm(A - Q*R, 'fro') / norm(A, 'fro');
R = randlapack.cqrrt(A); % R only: Q is not formed
```

CQRRT does not detect rank: on a rank-deficient input it usually returns without error, with a numerically singular R (it stops with an error only on an exactly singular sketch, e.g. a zero column, or a Cholesky failure). If the rank is uncertain, check `abs(diag(R))` or use `randlapack.cqrrpt`. Both take `'d_factor'` and `'state'` name-value options (CQRRPT also `'eps'`); see `help randlapack.cqrrpt` and `help randlapack.cqrrt`.

## Runtime notes

On **Linux MATLAB R2026a + Ubuntu 24.04** (and similar new-glibc hosts), MATLAB ships an older `libstdc++.so.6` than your compiler's. MEX files built against the system `libstdc++` fail at load with `version GLIBCXX_3.4.32 not found`. Workaround:
Expand All @@ -76,5 +99,5 @@ On **Windows**, the MEX uses MATLAB's own OpenMP runtime (`libiomp5md`) instead

## Conventions

* `J` is **1-based**, matching MATLAB. BQRRP stores 1-based pivots natively; no conversion in the MEX layer.
* `J` is **1-based**, matching MATLAB. BQRRP and CQRRPT store 1-based pivots natively; no conversion in the MEX layer.
* RNG `state` is Philox4x32: `state.counter` is `uint32[4]`, `state.key` is `uint32[2]`. Omit it (defaults to seed 0) or pass a scalar `uint32` seed.
36 changes: 3 additions & 33 deletions matlab/+randlapack/bqrrp.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@
%
% Two output modes, selected via a trailing string argument:
%
% 'explicit' (default): matches MATLAB qr(A, 'vector')
% 'explicit' (default): matches MATLAB qr(A, 'econ', 'vector')
% [Q, R, J] = randlapack.bqrrp(A)
% [Q, R, J] = randlapack.bqrrp(A, 'explicit')
% Q is m-by-k orthogonal, R is k-by-n upper-triangular, k = min(m, n),
Expand Down Expand Up @@ -50,7 +50,7 @@
% pass a smaller value to actually exercise the blocked structure.
b_sz = int64(min(m, n));
d_factor = cast(1.25, class(A));
state = make_default_state();
state = rl_default_state();
mode = 'explicit';

% --- Argument parsing ---
Expand All @@ -75,7 +75,7 @@
d_factor = cast(positional{2}, class(A));
end
if numel(positional) >= 3 && ~isempty(positional{3})
state = normalize_state(positional{3});
state = rl_normalize_state(positional{3}, mfilename);
end
if numel(positional) >= 4
error('randlapack:bqrrp:nargin', ...
Expand Down Expand Up @@ -104,33 +104,3 @@
% out2 = tau (length-n Householder scalars)
[out1, out2, J, state_out] = bqrrp_mex(A, b_sz, d_factor, state, mode);
end


function state = make_default_state()
state.counter = uint32([0, 0, 0, 0]);
state.key = uint32([0, 0]);
end


function state = normalize_state(in)
if isstruct(in)
if ~isfield(in, 'counter') || ~isfield(in, 'key')
error('randlapack:bqrrp:state', ...
'state struct must have fields .counter and .key');
end
validateattributes(in.counter, {'uint32'}, {'vector', 'numel', 4}, ...
mfilename, 'state.counter');
validateattributes(in.key, {'uint32'}, {'vector', 'numel', 2}, ...
mfilename, 'state.key');
state.counter = uint32(in.counter(:)).';
state.key = uint32(in.key(:)).';
elseif isnumeric(in) && isscalar(in)
% Treat scalar as a seed; place in key[1]. Counter starts at zero.
seed = uint32(in);
state.counter = uint32([0, 0, 0, 0]);
state.key = uint32([0, seed]);
else
error('randlapack:bqrrp:state', ...
'state must be a struct or a scalar seed');
end
end
124 changes: 124 additions & 0 deletions matlab/+randlapack/cqrrpt.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,124 @@
function [Q, R, J, state_out] = cqrrpt(A, varargin)
%RANDLAPACK.CQRRPT Randomized rank-revealing QR with column pivoting.
%
% [Q, R, J] = randlapack.cqrrpt(A) computes a rank-revealing QR
% factorization of a TALL matrix A, so that
% A(:, J) ~= Q * R
% with Q m-by-k orthonormal, R k-by-n upper-trapezoidal, and J a length-n
% vector of 1-based column pivots. Output layout matches
% qr(A, 'econ', 'vector').
%
% k is DISCOVERED, not requested: CQRRPT reveals the numerical rank of A
% and that rank is the inner dimension of Q and R. On a full-rank input
% k = n. Read it off size(R, 1).
%
% What the factorization guarantees, and what it does not: Q is
% orthonormal and the LEADING k pivoted columns always reproduce exactly,
% A(:, J(1:k)) == Q * R(:, 1:k)
% to machine precision. The full identity A(:, J) == Q*R additionally
% needs k to have reached the numerical rank of A: the ordinary case,
% including on rank-deficient input. It can fall short when CQRRPT's
% internal Cholesky QR breaks down and the conservative fallback estimate
% takes over, so check size(R, 1) and the residual on inputs whose
% conditioning you do not know.
%
% CQRRPT sketches A down to d = d_factor*n rows, so it requires
% m >= d_factor*n. For wide or square matrices use randlapack.bqrrp.
%
% Name-value options:
% 'eps' conditioning tolerance for rank detection. Default:
% sqrt(eps(class(A))).
%
% eps is NOT a singular-value truncation threshold.
% CQRRPT first truncates using a fixed machine-epsilon
% cutoff, so columns well above machine precision are
% kept whatever eps says. eps takes effect only in the
% re-estimation that runs when the internal Cholesky QR
% fails (e.g. on an ill-conditioned input), where the rank
% is cut once the R diagonal spans a ratio of
% sqrt(eps/eps(class(A))). It therefore bounds the
% orthogonality loss, which for Cholesky QR scales like
% u*cond(R)^2. Smaller eps means a stricter bound and a
% smaller returned rank.
% 'd_factor' sketch embedding factor (>= 1.0). Default: 1.25.
% Larger values sketch more rows: steadier rank
% detection, more work, and a taller input requirement.
% 'state' RNG state: a uint32 scalar seed, or a struct with
% .counter (uint32[4]) and .key (uint32[2]).
%
% [..., state] = randlapack.cqrrpt(...) returns the advanced RNG state.
%
% A must be single or double precision, real, and 2D.
%
% Example:
% B = randn(5000, 40);
% A = [B, B(:, 1:10)]; % 50 columns, rank 40
% [Q, R, J] = randlapack.cqrrpt(A);
% detected_rank = size(R, 1) % ~40
% err = norm(A(:, J) - Q*R, 'fro') / norm(A, 'fro');
%
% See also QR, RANDLAPACK.BQRRP, RANDLAPACK.CQRRT.

% --- Required argument check ---
validateattributes(A, {'single', 'double'}, {'2d', 'real', 'finite'}, ...
mfilename, 'A', 1);
[m, n] = size(A);
if m == 0 || n == 0
error('randlapack:cqrrpt:emptyInput', 'A must be nonempty');
end

% --- Defaults ---
% eps default: sqrt of machine epsilon in the working precision. eps
% bounds the orthogonality loss in the fallback rank estimate (see the
% help above), and Cholesky QR loses orthogonality like u*cond(R)^2, so
% sqrt(u) is its natural scale.
opts = struct('eps', sqrt(eps(class(A))), ...
'd_factor', 1.25, ...
'state', rl_default_state());

% --- Name-value parsing ---
if mod(numel(varargin), 2) ~= 0
error('randlapack:cqrrpt:nargin', ...
'Options must be given as name-value pairs.');
end
known = fieldnames(opts);
for i = 1:2:numel(varargin)
name = varargin{i};
if ~(ischar(name) || (isstring(name) && isscalar(name)))
error('randlapack:cqrrpt:optionName', ...
'Option names must be character vectors or strings.');
end
name = char(name);
idx = find(strcmpi(known, name), 1);
if isempty(idx)
error('randlapack:cqrrpt:unknownOption', ...
'Unknown option ''%s''. Valid options: %s.', ...
name, strjoin(known', ', '));
end
opts.(known{idx}) = varargin{i + 1};
end

% --- Validation ---
validateattributes(opts.eps, {'numeric'}, {'scalar', 'real', 'nonnegative'}, ...
mfilename, 'eps');
validateattributes(opts.d_factor, {'numeric'}, {'scalar', 'real', '>=', 1}, ...
mfilename, 'd_factor');

% The tall-matrix requirement is checked here as well as in the MEX, so
% the common mistake is caught before any marshalling happens and the
% message can name the MATLAB-level alternative.
d = floor(opts.d_factor * n);
if m < d
error('randlapack:cqrrpt:notTall', ...
['CQRRPT requires m >= d_factor*n. Got m=%d, n=%d, ' ...
'd_factor=%g (needs m >= %d). Lower d_factor toward 1.0, ' ...
'or use randlapack.bqrrp, which has no tall requirement.'], ...
m, n, opts.d_factor, d);
end

state = rl_normalize_state(opts.state, mfilename);

% --- Dispatch to MEX ---
[Q, R, J, state_out] = cqrrpt_mex(A, cast(opts.d_factor, class(A)), ...
cast(opts.eps, class(A)), state);
end
109 changes: 109 additions & 0 deletions matlab/+randlapack/cqrrt.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
function [out1, R, state_out] = cqrrt(A, varargin)
%RANDLAPACK.CQRRT Randomized Cholesky QR of a tall full-rank matrix.
%
% [Q, R] = randlapack.cqrrt(A) computes an unpivoted QR factorization of a
% TALL, FULL-RANK matrix A, so that
% A ~= Q * R
% with Q m-by-n orthonormal and R n-by-n upper-triangular. Output layout
% matches qr(A, 0) (the economy-size QR).
%
% R = randlapack.cqrrt(A), with a single output, returns R only and skips
% forming Q (CQRRT's R-only mode), like R = qr(A, 0). This is the cheaper
% call when only R is needed, e.g. to solve least-squares problems or to
% precondition with R.
%
% How it works: CQRRT sketches A down to d = d_factor*n rows, takes a QR of
% the small sketch, uses its R factor to precondition A, and finishes with
% Cholesky QR on the well-conditioned result. It does not pivot and does
% NOT detect rank: on a rank-deficient A it usually returns without error,
% with Q*R still reproducing A but R numerically singular (tiny diagonal
% entries). It stops with an error only when the sketch's R is exactly
% singular (for example, A has a zero column) or the Cholesky step fails.
% If the rank of A is uncertain, check abs(diag(R)) or use
% randlapack.cqrrpt, which reveals the rank (or randlapack.bqrrp).
%
% CQRRT requires m >= d_factor*n. For wide or square matrices use
% randlapack.bqrrp.
%
% Name-value options:
% 'd_factor' sketch embedding factor (>= 1.0). Default: 1.25.
% Larger values sketch more rows: a better
% preconditioner, more work, and a taller input
% requirement.
% 'state' RNG state: a uint32 scalar seed, or a struct with
% .counter (uint32[4]) and .key (uint32[2]).
%
% [Q, R, state] = randlapack.cqrrt(...) returns the advanced RNG state.
% (The single-output R-only form does not return it.)
%
% A must be single or double precision, real, and 2D.
%
% Example:
% A = randn(20000, 200);
% [Q, R] = randlapack.cqrrt(A);
% err = norm(A - Q*R, 'fro') / norm(A, 'fro');
% R_only = randlapack.cqrrt(A); % no Q formed
%
% See also QR, RANDLAPACK.CQRRPT, RANDLAPACK.BQRRP.

% --- Required argument check ---
validateattributes(A, {'single', 'double'}, {'2d', 'real', 'finite'}, ...
mfilename, 'A', 1);
[m, n] = size(A);
if m == 0 || n == 0
error('randlapack:cqrrt:emptyInput', 'A must be nonempty');
end

% --- Defaults ---
opts = struct('d_factor', 1.25, ...
'state', rl_default_state());

% --- Name-value parsing ---
if mod(numel(varargin), 2) ~= 0
error('randlapack:cqrrt:nargin', ...
'Options must be given as name-value pairs.');
end
known = fieldnames(opts);
for i = 1:2:numel(varargin)
name = varargin{i};
if ~(ischar(name) || (isstring(name) && isscalar(name)))
error('randlapack:cqrrt:optionName', ...
'Option names must be character vectors or strings.');
end
name = char(name);
idx = find(strcmpi(known, name), 1);
if isempty(idx)
error('randlapack:cqrrt:unknownOption', ...
'Unknown option ''%s''. Valid options: %s.', ...
name, strjoin(known', ', '));
end
opts.(known{idx}) = varargin{i + 1};
end

% --- Validation ---
validateattributes(opts.d_factor, {'numeric'}, {'scalar', 'real', '>=', 1}, ...
mfilename, 'd_factor');

% The tall-matrix requirement is checked here as well as in the MEX, so
% the common mistake is caught before any marshalling happens and the
% message can name the MATLAB-level alternative.
d = floor(opts.d_factor * n);
if m < d
error('randlapack:cqrrt:notTall', ...
['CQRRT requires m >= d_factor*n. Got m=%d, n=%d, ' ...
'd_factor=%g (needs m >= %d). Lower d_factor toward 1.0, ' ...
'or use randlapack.bqrrp, which has no tall requirement.'], ...
m, n, opts.d_factor, d);
end

state = rl_normalize_state(opts.state, mfilename);

% --- Dispatch to MEX ---
d_factor = cast(opts.d_factor, class(A));
if nargout <= 1
% R-only mode: Q is not formed.
[~, out1] = cqrrt_mex(A, d_factor, false, state);
else
[out1, R, state_out] = cqrrt_mex(A, d_factor, true, state);
end
end
11 changes: 11 additions & 0 deletions matlab/+randlapack/private/rl_default_state.m
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,11 @@
function state = rl_default_state()
%RL_DEFAULT_STATE Default Philox4x32 RNG state (all zeros).
%
% Shared by the randlapack.* wrappers. Functions in a private/ directory
% are visible to the enclosing package's functions and to nobody else,
% which is exactly the scope these helpers want. The compiled MEX files
% land in this same directory, for the same reason.

state.counter = uint32([0, 0, 0, 0]);
state.key = uint32([0, 0]);
end
Loading