What does the lu function do in MATLAB / RunMat?

lu(A) computes the LU factorization of a real or complex matrix A using partial pivoting. It exposes the same calling forms as MATLAB:

• Single output: lu(A) returns a single matrix whose strictly lower-triangular entries encode L (with an implicit unit diagonal) and whose upper-triangular part encodes U.
• Two outputs: [L, U] = lu(A) returns the explicit unit-lower-triangular factor L and the upper-triangular factor U.
• Three outputs: [L, U, P] = lu(A) additionally returns a permutation so that P * A = L * U. Use the option 'vector' to receive the permutation as a pivot vector instead of a matrix.

The implementation follows MATLAB’s dense behaviour for full matrices and supports rectangular inputs.

How does the lu function behave in MATLAB / RunMat?

• Partial pivoting is applied to improve numerical stability. The permutation is encoded either as a dense matrix ('matrix', default) or as a pivot vector ('vector').
• Rectangular inputs are supported. L is always m × m (unit lower-triangular), and U is m × n, where m and n are the row and column counts of A.
• Singular matrices are permitted. Zero pivots propagate into the U factor just as in MATLAB; MATLAB-compatible warnings are not yet emitted.
• Only the first three outputs are implemented today. Column permutations (Q) and scaling (R) for the five-output sparse form are not yet available.

GPU execution in RunMat

• When an acceleration provider implements the lu hook (the WGPU provider does), the factorization executes through that provider and the combined LU factor, L, U, and permutation outputs all remain on the device. The current WGPU backend performs the decomposition on the host once and immediately reuploads the factors so residency is preserved until dedicated kernels land.
• The 'vector' option likewise returns a GPU-resident pivot vector when a provider hook is active.
• If no provider hook is available, RunMat automatically gathers the input to host memory and falls back to the CPU implementation so behaviour stays MATLAB-compatible.

Examples of using the lu function in MATLAB / RunMat

Factorizing a square matrix with lu

A = [2 1 1; 4 -6 0; -2 7 2];
[L, U, P] = lu(A);

Expected output (up to floating-point roundoff):

L =
     1     0     0
    -1     1     0
     0    -1     1

U =
     4    -6     0
     0     1     1
     0     0     3

P =
     0     1     0
     1     0     0
     0     0     1

Obtaining only the combined LU factor

LU = lu([1 3 5; 2 4 7; 1 1 0]);

Expected output:

LU =
     2     4     7
   0.5     1   -1.5
   0.5   -0.5    2

Requesting the permutation vector with the 'vector' option

[L, U, p] = lu([4 3; 6 3], 'vector');

Expected output:

p =
     2
     1

LU factorization of a rectangular matrix

A = [3 1 2; 6 3 4];
[L, U, P] = lu(A);

Expected output:

L =
     1     0
     0.5   1

U =
     6     3     4
     0    -0.5    0

P =
     0     1
     1     0

Using LU factors to solve a linear system

A = [3 1 2; 6 3 4];
b = [1; 2];
[L, U, P] = lu(A);
y = L \ (P * b);
x = U \ y;

Expected output:

x =
    0.0
    0.5
   -0.0

Running lu on a gpuArray

G = gpuArray([10 7; 3 2]);
[L, U, P] = lu(G);
class(L)
class(U)
class(P)

Expected output:

ans =
    'gpuArray'

ans =
    'gpuArray'

ans =
    'gpuArray'

If no acceleration provider exposes lu, RunMat gathers the input and returns the factors as host double arrays instead.

FAQ

Why does RunMat currently stop at three outputs?

Column pivoting (Q) and scaling (R) from MATLAB’s five-output sparse form are planned but not yet implemented. The dense three-output contract mirrors MATLAB’s default dense behaviour.

Does the permutation vector use MATLAB’s 1-based indexing?

Yes. When you request 'vector', the returned pivot vector contains 1-based row indices so that A(p, :) = L * U.

How are singular matrices handled?

Partial pivoting proceeds exactly as in MATLAB. If a pivot column is entirely zero, the corresponding diagonal entries in U become zero. No warning is emitted yet.

Are complex matrices supported?

Yes. Complex inputs produce complex L, U, and LU. The permutation remains real because it only contains zeros and ones.

Will the factors stay on the GPU when I pass a gpuArray?

Yes. When the active acceleration provider exposes the lu hook (WGPU today), the combined factor, L, U, and the permutation outputs remain gpuArray values—the provider currently performs the decomposition on the host once and reuploads the results to preserve residency. Without provider support, RunMat gathers to host memory before returning the factors.

Can I call lu on logical arrays?

Yes. Logical inputs are promoted to double precision before factorization, matching MATLAB semantics.

Is pivoting deterministic?

Yes. Partial pivoting always chooses the first maximal entry in each column, mirroring MATLAB’s behaviour for dense matrices.

How accurate is the factorization?

The implementation uses standard double-precision arithmetic (or complex double when needed). Numerical properties therefore match MATLAB’s dense fallback (without iterative refinement).

What happens if I pass more than one option argument?

RunMat currently supports at most one option string ('matrix' or 'vector'). Passing additional options raises an error.

Can I reuse the combined LU factor to solve systems?

Yes. The combined matrix returned by lu(A) stores L in the strictly lower-triangular part (with an implicit unit diagonal) and U in the upper-triangular part, just like MATLAB. You can use forward/back substitution routines that understand this layout.

See Also

det, inv, chol, qr, solve, gpuArray

Source & Feedback

• Implementation: crates/runmat-runtime/src/builtins/math/linalg/factor/lu.rs
• Found an issue or missing behaviour? Open an issue with details and a minimal reproduction.