API Reference

nmf(X; rank::Int=10, maxiter::Int=500, tol::Real=1e-4, seed=nothing)

Compute a standard non-negative matrix factorization (NMF) of X using multiplicative update rules with squared Frobenius reconstruction loss.

Arguments

• X::AbstractMatrix{<:Real}: Non-negative data matrix of size (m, n).

Keyword Arguments

• rank::Int=10: Target factorization rank.
• maxiter::Int=500: Maximum number of iterations.
• tol::Real=1e-4: Relative tolerance for stopping based on objective change.
• seed=nothing: Optional random seed for reproducible initialization.

Returns

• W::Matrix{Float64}: Non-negative basis matrix of size (m, rank).
• H::Matrix{Float64}: Non-negative coefficient matrix of size (rank, n).
• history::Vector{Float64}: Squared Frobenius reconstruction error per iteration.

Side Effects

• None. (The function does not modify X.)

Errors

• None.

Notes

• Uses multiplicative updates with a small epsilon for numerical stability.
• Convergence check is based on relative change in the objective.
• The reconstruction can be obtained as W * H.
• Input data should be non-negative.

Examples

julia> X, _, _ = generate_synthetic_data(20, 15; rank=3, seed=42);

julia> W, H, history = nmf(X; rank=3, maxiter=200, tol=1e-5, seed=42);

julia> size(W), size(H), length(history) > 0
((20, 3), (3, 15), true)

robustnmf(X; kwargs...)

Default robust NMF entry point.

This calls Huber-loss robust NMF implementation by default.

Arguments

• X::AbstractMatrix{<:Real}: Non-negative data matrix (m, n).

Keyword Arguments

• kwargs...: Forwarded to robustnmf_huber.

Returns

• (W, H, history): See robustnmf_huber.

Side Effects

• None.

Errors

• ArgumentError: If inputs are invalid (see robustnmf_huber).

Notes

• Public API entry point; keeps the external name stable if the default robust method changes.

Examples

julia> using RobustNMF

julia> X, _, _ = generate_synthetic_data(20, 12; rank=4, seed=7);

julia> W, H, history = robustnmf(X; rank=4, maxiter=50, tol=1e-6, delta=1.0, seed=7);

julia> size(W), size(H)
((20, 4), (4, 12))

robustnmf_huber(X; rank=10, maxiter=500, tol=1e-4, delta=1.0, seed=nothing)

Implements a robust NMF variant using the Huber loss on the reconstruction residuals: R = X - W*H

Huber loss (element-wise) with threshold δ:

• if |r| ≤ δ: 0.5 * r^2
• else: δ * (|r| - 0.5*δ)

We optimize it using an IRLS-style weighted least squares approach:

1. compute residual R
2. compute weights Ω = huber_weights(R, δ)
3. perform weighted multiplicative updates (update_huber)
4. track huber_loss(R, δ) in history
5. stop when relative change in objective is below tol

Arguments

• X::AbstractMatrix{<:Real}: Non-negative data matrix of size (m, n).

Keyword Arguments

• rank::Int=10: Factorization rank.
• maxiter::Int=500: Maximum number of iterations.
• tol::Real=1e-4: Relative tolerance for stopping based on objective change.
• delta::Real=1.0: Huber threshold δ (must be > 0).
• seed=nothing: Optional random seed for reproducibility.

Returns

• W::Matrix{Float64}: Non-negative basis matrix of size (m, rank).
• H::Matrix{Float64}: Non-negative coefficient matrix of size (rank, n).
• history::Vector{Float64}: Huber objective values per iteration.

Side Effects

• None. (The function does not modify X.)

Errors

• ArgumentError: If X contains negative entries or parameters are invalid.

Notes

• Uses IRLS weights via huber_weights and weighted multiplicative updates.
• Convergence check is based on relative change of Huber objective.
• Initialization uses a local RNG seeded by seed when provided.

Examples

julia> using RobustNMF

julia> X, _, _ = generate_synthetic_data(30, 20; rank=5, seed=123);

julia> W, H, history = robustnmf_huber(X; rank=5, maxiter=100, tol=1e-6, delta=1.0, seed=123);

julia> size(W), size(H), length(history) > 0
((30, 5), (5, 20), true)

robustnmf_l21(X; rank=10, maxiter=500, tol=1e-4, seed=nothing)

L2,1-Norm Regularized Non-negative Matrix Factorization.

Minimizes: ||X - FG||_{2,1} where the L2,1-norm promotes robustness to sample-wise outliers (entire corrupted columns in X).

Arguments

• X::AbstractMatrix{<:Real}: Non-negative data matrix (m, n).

Keyword Arguments

• rank::Int=10: Factorization rank.
• maxiter::Int=500: Maximum number of iterations.
• tol::Real=1e-4: Absolute tolerance for stopping.
• seed=nothing: Optional random seed for reproducibility.

Returns

• F::Matrix{Float64}: Non-negative basis matrix (m, rank).
• G::Matrix{Float64}: Non-negative coefficient matrix (rank, n).
• history::Vector{Float64}: L2,1 objective values per iteration.

Side Effects

• None. (The function does not modify X.)

Errors

• AssertionError: If X contains negative entries or parameters are invalid.

Notes

• Legacy algorithm kept for compatibility; Huber is the default robust method.
• Convergence check uses absolute objective value (error < tol).

Examples

julia> using RobustNMF

julia> X, _, _ = generate_synthetic_data(20, 12; rank=4, seed=7);

julia> F, G, history = robustnmf_l21(X; rank=4, maxiter=50, tol=1e-3, seed=7);

julia> size(F), size(G), length(history) > 0
((20, 4), (4, 12), true)

update_huber(X, W, H, Ω; ϵ=eps(Float64))

Perform one weighted multiplicative update step for robust NMF with Huber IRLS weights.

We interpret Huber as a weighted least-squares problem at each IRLS step:

min{W,H ≥ 0} ||Ω ⊙ (X - W*H)||_F^2

Given Ω (same size as X), the standard Frobenius multiplicative updates become: H ← H ⊙ (W' * (Ω ⊙ X)) ./ (W' * (Ω ⊙ (WH)) + ϵ) W ← W ⊙ ((Ω ⊙ X) * H') ./ ((Ω ⊙ (WH)) * H' + ϵ)

Returns updated (W, H).

Arguments

• X::AbstractMatrix{<:Real}: Non-negative data matrix (m, n).
• W::AbstractMatrix{<:Real}: Current basis matrix (m, rank).
• H::AbstractMatrix{<:Real}: Current coefficient matrix (rank, n).
• Ω::AbstractMatrix{<:Real}: IRLS weight matrix (same size as X).

Keyword Arguments

• ϵ=eps(Float64): Small constant for numerical stability.

Returns

• W::AbstractMatrix{<:Real}: Updated basis matrix.
• H::AbstractMatrix{<:Real}: Updated coefficient matrix.

Side Effects

• Updates W and H in-place.

Errors

• None.

Notes

• Multiplicative updates preserve non-negativity if W and H start non-negative.
• Uses Ω ⊙ X and Ω ⊙ (W*H) to compute weighted updates.

huber_loss(R::AbstractMatrix{<:Real}, delta::Real; ϵ=eps(Float64))

Compute the Huber loss for a residual matrix R.

Huber loss is robust to outliers:

• For small residuals it behaves like squared error (L2).
• For large residuals it behaves like absolute error (L1), reducing the impact of outliers.

Element-wise definition for residual r:

• if |r| ≤ δ: 0.5 * r^2
• else: δ * (|r| - 0.5*δ)

Returns the sum over all entries of R.

Arguments

• R::AbstractMatrix{<:Real}: Residual matrix.
• delta::Real: Huber threshold δ (> 0).

Keyword Arguments

• ϵ=eps(Float64): Small constant for numerical stability (used in validation/consistency).

Returns

• loss::Float64: Sum of Huber losses over all entries of R.

Side Effects

• None.

Errors

• ArgumentError: If delta <= 0.

Notes

• For |r| ≤ δ, the loss is quadratic; for |r| > δ, it is linear.
• Summation uses the element type of R for type stability.

huber_weights(R::AbstractMatrix{<:Real}, delta::Real; ϵ=eps(Float64))

Compute the Huber IRLS weights matrix Ω for a residual matrix R.

We use an iteratively reweighted least squares (IRLS) interpretation:

• Small residuals get weight 1.0 (quadratic region).
• Large residuals get weight δ / (|r| + ϵ), which downweights outliers.

Element-wise:

• if |r| ≤ δ: w = 1
• else: w = δ / (|r| + ϵ)

Returns Ω with the same size as R.

Arguments

• R::AbstractMatrix{<:Real}: Residual matrix.
• delta::Real: Huber threshold δ (> 0).

Keyword Arguments

• ϵ=eps(Float64): Small constant to avoid division by zero.

Returns

• Ω::Matrix{Float64}: Weight matrix of same size as R.

Side Effects

• None.

Errors

• ArgumentError: If delta <= 0.

Notes

• Weights are 1.0 in the quadratic region and decrease for large residuals.
• Used to build weighted multiplicative updates in IRLS.

update_l21(X, F, G; eps_update=1e-10)

Perform one iteration of L2,1-NMF multiplicative updates.

The L2,1-norm promotes row sparsity in the residual matrix, making the algorithm robust to sample-wise (column-wise) outliers.

Arguments

• X::AbstractMatrix: Data matrix (m, n).
• F::AbstractMatrix: Current basis matrix (m, rank).
• G::AbstractMatrix: Current coefficient matrix (rank, n).

Keyword Arguments

• eps_update::Real=1e-10: Small constant for numerical stability.

Returns

• F_new::Matrix{Float64}: Updated basis matrix.
• G_new::Matrix{Float64}: Updated coefficient matrix.

Side Effects

• None.

Errors

• None.

Notes

• Legacy algorithm kept for compatibility; Huber is the default robust method.
• Update uses a diagonal reweighting matrix D derived from column residuals.

l21_loss(X::AbstractMatrix)

Compute the L2,1-norm of matrix X. The L2,1-norm is the sum of the L2-norms of each column.

Arguments

• X::AbstractMatrix: Input matrix.

Keyword Arguments

• None.

Returns

• value::Float64: Sum of L2-norms of columns.

Side Effects

• None.

Errors

• None.

Notes

• Legacy helper used by L2,1-NMF.

Examples

julia> X = [1.0 2.0; 3.0 4.0];

julia> l21_loss(X) >= 0
true

generate_synthetic_data(m::Int, n::Int; rank::Int=10, noise_level::Real=0.0, seed=nothing)

Generate a synthetic non-negative data matrix X as W * H with random non-negative factors. Optionally adds Gaussian noise and clips negative values to keep X ≥ 0.

Arguments

• m::Int: Number of rows of X.
• n::Int: Number of columns of X.

Keyword Arguments

• rank::Int=10: Rank of the factorization.
• noise_level::Real=0.0: Standard deviation of Gaussian noise.
• seed=nothing: Optional random seed for reproducibility.

Returns

• X::Matrix{Float64}: Generated non-negative data matrix.
• W::Matrix{Float64}: Left factor of size (m, rank).
• H::Matrix{Float64}: Right factor of size (rank, n).

Side Effects

• None.

Errors

• None.

Notes

• Uses a local RNG seeded with seed (if provided) for deterministic output.
• When noise_level > 0, the result is clipped at 0.0 to enforce non-negativity.

Examples

julia> X, W, H = generate_synthetic_data(20, 15; rank=3, seed=42);

julia> size(W), size(H)
((20, 3), (3, 15))

julia> minimum(X) >= 0
true

add_gaussian_noise!(X::AbstractMatrix; σ::Real=0.1, clip_at_zero::Bool=true)

Add Gaussian noise with standard deviation σ to the matrix X in-place. Optionally clip negative entries to preserve non-negativity.

Arguments

• X::AbstractMatrix: Data matrix to be corrupted.

Keyword Arguments

• σ::Real=0.1: Noise standard deviation.
• clip_at_zero::Bool=true: Enforce non-negativity after corruption.

Returns

• X: The modified input matrix (in-place).

Side Effects

• Modifies X in-place.

Errors

• None.

Notes

• Uses randn! to generate Gaussian noise.
• If clip_at_zero=true, negative values are replaced by 0.0.

Examples

julia> X = abs.(randn(5, 5));

julia> add_gaussian_noise!(X; σ=0.2);

julia> minimum(X) >= 0
true

add_sparse_outliers!(X::AbstractMatrix; fraction::Real=0.01, magnitude::Real=5.0, seed=nothing)

Add sparse, large positive outliers to a fraction of the entries of X in-place.

Arguments

• X::AbstractMatrix: Data matrix to be corrupted.

Keyword Arguments

• fraction::Real=0.01: Fraction of entries to corrupt.
• magnitude::Real=5.0: Maximum outlier amplitude.
• seed=nothing: Optional random seed.

Returns

• X: The modified input matrix (in-place).

Side Effects

• Modifies X in-place.

Errors

• None.

Notes

• Corrupts max(1, round(Int, fraction * m*n)) entries.
• Added outliers are sampled from Uniform(0, magnitude).

Examples

julia> X = zeros(10, 10);

julia> add_sparse_outliers!(X; fraction=0.05, seed=1);

julia> count(>(0.0), X) > 0
true

normalize_nonnegative!(X::AbstractMatrix; rescale::Bool=true)

Shift X in-place so that the minimum becomes 0.0, and optionally rescale to [0, 1].

Arguments

• X::AbstractMatrix: Input matrix.

Keyword Arguments

• rescale::Bool=true: Whether to divide by the maximum value after shifting.

Returns

• X: The normalized matrix (in-place).

Side Effects

• Modifies X in-place.

Errors

• None.

Notes

• If rescale=true and maximum(X) == 0, rescaling is skipped.

Examples

julia> X = [-1.0 2.0; 3.0 -4.0];

julia> normalize_nonnegative!(X);

julia> minimum(X), maximum(X)
(0.0, 1.0)

load_image_folder(dir::AbstractString; pattern::AbstractString=".png", normalize::Bool=true)

Load images from a folder, convert to grayscale, flatten, and stack them as columns of X.

Arguments

• dir::AbstractString: Path to the image directory.

Keyword Arguments

• pattern::AbstractString=".png": File extension filter (matched via endswith).
• normalize::Bool=true: Normalize output matrix to [0, 1].

Returns

• X::Matrix{Float64}: One column per image.
• (height, width): Original image dimensions.
• filenames::Vector{String}: Loaded base file names.

Side Effects

• Reads image files from disk.

Errors

• ErrorException: If the directory does not exist or no files match pattern.
• ErrorException: If images have inconsistent sizes.

Notes

• Images are converted to grayscale and stored as Float64.
• If normalize=true, normalize_nonnegative! is applied to X.

Examples

julia> # X, size, names = load_image_folder("faces/")

plot_convergence(history::Vector;
                 title::String="NMF Convergence",
                 objective::Symbol=:auto,
                 ylabel::Union{Nothing,String}=nothing,
                 log_scale::Bool=true)

Plot the convergence history (history) produced by an NMF routine.

This function is intentionally objective-agnostic: depending on the algorithm, history may contain Frobenius reconstruction error (standard NMF), Huber loss (robust NMF), or another objective value.

Arguments

• history::Vector: Objective values recorded per iteration.

Keyword Arguments

• title::String: Plot title.
• objective::Symbol=:auto: Hint for labeling the objective.
  • :frobenius → "Frobenius Error"
  • :huber → "Huber Loss"
  • :l21 → "L2,1 Loss"
  • :auto → "Objective" (neutral default)
• ylabel::Union{Nothing,String}=nothing: Explicit y-axis label. If provided, this overrides objective.
• log_scale::Bool=true: Use logarithmic scale for y-axis.

Returns

• p::Plots.Plot: Plot object.

Side Effects

• None.

Errors

• None.

Notes

• Log scaling helps visualize early convergence behavior.

Examples

julia> using RobustNMF, Plots

julia> X, _, _ = generate_synthetic_data(20, 12; rank=4, seed=4);

julia> _, _, history = nmf(X; rank=4, maxiter=500, tol=1e-5, seed=101);

julia> p = plot_convergence(history; objective=:frobenius);

julia> _, _, history = robustnmf(X; rank=10, maxiter=500);

julia> plot_convergence(history; objective=:huber);

julia> p isa Plots.Plot
true

plot_basis_vectors(W::AbstractMatrix; img_shape=nothing, max_components::Int=16,
                   title::String="Basis Vectors (W)", layout=nothing)

Visualize the basis vectors (columns of W) as heatmaps or images.

Arguments

• W::AbstractMatrix: Basis matrix of size (m, rank).

Keyword Arguments

• img_shape: Tuple (height, width) to reshape each basis vector as an image. If nothing, displays as 1D heatmaps.
• max_components::Int=16: Maximum number of components to display.
• title::String: Plot title.
• layout: Custom layout tuple (rows, cols). If nothing, auto-computed.

Returns

• p::Plots.Plot: Plot object showing the basis vectors.

Side Effects

• None.

Errors

• ErrorException: If img_shape does not match the length of a basis vector.

Notes

• The number of displayed components is min(rank, max_components).
• Layout is chosen to be roughly square when not provided.

Examples

julia> using RobustNMF, Plots

julia> X, _, _ = generate_synthetic_data(40, 30; rank=6, seed=1);

julia> W, _, _ = nmf(X; rank=6, maxiter=50, tol=1e-5, seed=42);

julia> p = plot_basis_vectors(W; max_components=6);

julia> p isa Plots.Plot
true

plot_reconstruction_comparison(X_original::AbstractMatrix, X_recon::AbstractMatrix;
                               img_shape=nothing, n_samples::Int=5,
                               title::String="Reconstruction Comparison")

Compare original data with reconstructed data side by side.

Arguments

• X_original::AbstractMatrix: Original data matrix.
• X_recon::AbstractMatrix: Reconstructed data matrix (W * H).

Keyword Arguments

• img_shape: Tuple (height, width) for reshaping columns as images.
• n_samples::Int=5: Number of samples to display.
• title::String: Plot title.

Returns

• p::Plots.Plot: Plot object showing original vs reconstructed samples.

Side Effects

• None.

Errors

• AssertionError: If X_original and X_recon have different dimensions.
• ErrorException: If img_shape is incompatible with column length.

Notes

• For image data, each sample is displayed as an image pair.
• For vector data, each sample is plotted as a line comparison.

Examples

julia> using RobustNMF, Plots

julia> X, _, _ = generate_synthetic_data(20, 12; rank=4, seed=3);

julia> W, H, _ = nmf(X; rank=4, maxiter=50, tol=1e-5, seed=100);

julia> p = plot_reconstruction_comparison(X, W * H; n_samples=3);

julia> p isa Plots.Plot
true

plot_nmf_summary(X::AbstractMatrix,
                 W::AbstractMatrix,
                 H::AbstractMatrix,
                 history::Vector;
                 img_shape=nothing,
                 max_basis::Int=9,
                 max_samples::Int=4,
                 objective::Symbol=:auto,
                 convergence_ylabel::Union{Nothing,String}=nothing,
                 title::String="NMF Summary")

Create a summary visualization for an NMF result.

The summary consists of:

1. Basis vectors / components (columns of W)
2. Activating coefficients (rows of H)
3. Reconstruction comparison (original vs reconstructed data)
4. Convergence curve (objective value over iterations)

This function is algorithm-agnostic and works for:

• Standard NMF (Frobenius objective)
• Robust NMF with Huber loss
• Legacy L2,1 NMF

Arguments

• X::AbstractMatrix: Original non-negative data matrix (m × n).
• W::AbstractMatrix: Basis matrix (m × r).
• H::AbstractMatrix: Coefficient matrix (r × n).
• history::Vector: Objective values recorded during optimization.

Keyword Arguments

• img_shape: Optional tuple (height, width) if columns of X or W represent vectorized images.
• max_basis::Int=9: Maximum number of basis vectors (columns of W) to visualize.
• max_samples::Int=4: Maximum number of samples / activations to visualize.
• objective::Symbol=:auto: Type of objective used to generate history.
  • :frobenius → squared Frobenius reconstruction objective
  • :huber → Huber loss (robust NMF)
  • :l21 → L2,1 loss
  • :auto → neutral label ("Objective")
• convergence_ylabel::Union{Nothing,String}=nothing: Explicit y-axis label for the convergence plot. If provided, this overrides the label implied by objective.
• title::String="NMF Summary": Overall title for the summary figure.

Returns

• p::Plots.Plot: Plot object with a comprehensive summary.

Side Effects

• None.

Errors

• None.

Notes

• Includes a text panel with error metrics and iteration count.
• For image data, set img_shape to visualize basis vectors and reconstructions.

Examples

This example runs the summary once for standard NMF and once for robust NMF (Huber).

julia> using RobustNMF, Plots

julia> X, _, _ = generate_synthetic_data(30, 20; rank=5, seed=5);

julia> W, H, history = nmf(X; rank=5, maxiter=40, tol=1e-5, seed=102);

julia> p = plot_nmf_summary(X, W, H, history; objective=:frobenius, max_basis=4, max_samples=2);

julia> X, _, _ = generate_synthetic_data(30, 20; rank=5, seed=5);

julia> W, H, history = robustnmf(X; rank=5, maxiter=40, tol=1e-5, seed=103);

julia> p = plot_nmf_summary(X, W, H, history; objective=:huber, max_basis=4, max_samples=2);

julia> p isa Plots.Plot
true

plot_activation_coefficients(H::AbstractMatrix; max_samples::Int=10,
                             title::String="Activation Coefficients (H)")

Visualize the activation coefficient matrix H as a heatmap or as individual sample profiles.

Arguments

• H::AbstractMatrix: Coefficient matrix of size (rank, n).

Keyword Arguments

• max_samples::Int=10: Maximum number of samples to display (if showing individual profiles).
• title::String: Plot title.

Returns

• p::Plots.Plot: Plot object.

Side Effects

• None.

Errors

• None.

Notes

• If n ≤ 100 and rank ≤ 50, a full heatmap is shown.
• Otherwise, bar plots for up to max_samples samples are shown.

Examples

julia> using RobustNMF, Plots

julia> X, _, _ = generate_synthetic_data(30, 20; rank=5, seed=2);

julia> _, H, _ = nmf(X; rank=5, maxiter=50, tol=1e-5, seed=99);

julia> p = plot_activation_coefficients(H; max_samples=5);

julia> p isa Plots.Plot
true

plot_image_reconstruction(X::AbstractMatrix, W::AbstractMatrix, H::AbstractMatrix,
                          img_shape::Tuple{Int,Int}; indices=nothing, n_images::Int=5)

Specialized function for visualizing image reconstruction quality. Shows original, reconstructed, and difference images side by side.

Arguments

• X::AbstractMatrix: Original image data (each column is a flattened image).
• W::AbstractMatrix: Basis matrix.
• H::AbstractMatrix: Coefficient matrix.
• img_shape::Tuple{Int,Int}: Image dimensions (height, width).

Keyword Arguments

• indices: Specific image indices to display. If nothing, randomly selected.
• n_images::Int=5: Number of images to display.

Returns

• p::Plots.Plot: Plot object.

Side Effects

• None.

Errors

• ErrorException: If img_shape does not match column length.

Notes

• The difference image uses absolute error |X - W*H|.
• If indices is provided, only the first n_images indices are used.

Examples

julia> using RobustNMF, Plots

julia> X, _, _ = generate_synthetic_data(100, 8; rank=5, seed=6);

julia> W, H, _ = nmf(X; rank=5, maxiter=40, tol=1e-5, seed=104);

julia> p = plot_image_reconstruction(X, W, H, (10, 10); n_images=3);

julia> p isa Plots.Plot
true