Code

I have (finally) started to put my codes on gitlab; see https://gitlab.com/ngillis/ 

Non-symmetric semidefinite Procrustes problem

This Matlab code contains several algorithms for solving the non-symmetric positive semidefinite (PSD) Procrustes problem: given rectangular matrices X and B, find the matrix A that minimizes the Frobenius norm of AX-B and such that A+AT is PSD. It contains a new semi-analytical approach that is guaranteed to produce an optimal minimum-norm solution.

See N. Gillis and S. Sicilia, "Minimum-norm solutions of the non-symmetric semidefinite Procrustes problem", June 2024. [arXiv]

Target Stationary Distribution Problem

This Matlab code allows you to solve the Target Stationary Distribution Problem (TSDP): Given an n-by-n stochastic and irreducible matrix G, and a target distribution μ, find a correction matrix Δ such that G+Δ is stochastic and has μ as its stationary distribution, that is, (G+Δ)Tμ=μ, while minimizing the component-wise ell-1 norm of Δ. It is also possible to add a constraint on the support of Δ. In summary, this codes allows you to solve the following optimization problem: minΔ ||Δ||1 = ∑i,j |Δ(i,j)| such that G+Δ≥0, (G+Δ)Tμ=μ, Δ1n​=0, supp(Δ)⊂Ω for some given support Ω. 

See N. Gillis and P. Van Dooren, "Assigning Stationary Distributions to Sparse Stochastic Matrices", December 2023. [arXiv]

Boolean Matrix Factorization

This Matlab code provides the implementation of several algorithms, including two new effective ones, to solve the following problem: Given a binary matrix X(mxn) and a factorization rank r, compute the binary matrices W(mxr) and H(rxn) such that  ||X - min(WH,1)||F2  is minimized. 

See C. Kolomvakis, A. Vandaele and N. Gillis, "Algorithms for Boolean Matrix Factorization using Integer Programming", May 2023. [arXiv]

Nonlinear Matrix Decomposition with the ReLU function

This Matlab code provides the implementation of several algorithms, including two new effective ones, to solve the following problem: Given a (sparse) nonnegative matrix X(mxn) and a factorization rank r, compute W(mxr) and H(rxn) such that  ||X - max(WH,0)||F2  is minimized. 

See G. Seraghiti, A. Awari, A. Vandaele, M. Porcelli, N. Gillis, "Accelerated Algorithms for Nonlinear Matrix Decomposition with the ReLU function", May 2023. [arXiv]

Approximating the nearest Ω-Stable Matrix, Ω being an LMI region

This Matlab code can be used to compute the nearest Ω-stable matrix from an unstable one, that is, to compute a matrix whose eigenvalues are within Ω and that is as close as possible to the given input matrix. The set Ω, within the complex plane, is an LMI region parametrized using a dissipative-Hamiltonian (DH) characterization.  

See N. Choudhary, N. Gillis, P. Sharma, "Characterizing matrices with eigenvalues in an LMI region: A dissipative-Hamiltonian approach", October 2022. [arXiv]

Subspace Clustering Using Data Augmentation

This Matlab code  solves the subspace clustering problem using a new model that extends sparse subspace clustering (SSC) using data augmentation, in the unsupervised and semi-supervised settings; see 

M. Abdolali and N. Gillis, "Revisiting data augmentation for subspace clustering", July 2021. [arXiv]

Note: The code from the (short) ICASSP paper, "Subspace Clustering Using Unsupervised Data Augmentation", ICASSP 2022, 22-27 May 2022, Singapore, [doi], is available [here]. 

Deep MFs

This Matlab code allows you to compute deep matrix factorizations (MFs), decomposing a matrix X as WLHL...H2H1, with different constraints (sparsity, nonnegativity) and regularizers (minimum-volume); see 

P. De Handschutter and N. Gillis, "A consistent and flexible framework for deep matrix factorizations",  Pattern Recognition 134, 109102, 2023. [doi] [arXiv]

Partial identifiability for NMF

This Matlab code allows to check the partial identifiability conditions provided in the paper below. In a nutshell, given an Exact NMF X = W H,  it provides a subset of columns of W and of rows of H that are guaranteed to be unique (that is, identified in any other Exact NMF of X); see

N. Gillis and R. Rajkó, "Partial Identifiability for Nonnegative Matrix Factorization", Accepted in SIAM J. on Matrix Analysis and Applications, 2022. [arXiv] 

Smoothed separable NMF

This Matlab code can be used to tackle smoothed separable nonnegative matrix factorization (NMF): Given a nonnegative data matrix X and an integer r, find r index sets Ki containing p indices (i=1,2,...r) such that there exists a nonnegative matrix H with X ≈ mean(X(:,Ki)) H where mean(X(:,Ki)) is the average of the column of X(:,Ki); see 

N. Nadisic, N. Gillis and C. Kervazo, "Smoothed Separable Nonnegative Matrix Factorization", October 2021. [arXiv] 

Bounded Simplex-Structured Matrix Factorization

This Matlab code in [Matlab] and [Julia] solves Bounded Simplex-Structured Matrix Factorization (BSSMF): given an m-by-n matrix X, a factorization r, and an interval [a,b] in Rm,  it looks for an m-by-r matrix W and an r-by-n matrix H such that X ≈ WH where W(:,k) in [a,b] for all k, and H is column stochastic ; see the paper 

O. Vu Thanh, N. Gillis and F. Lecron, "Bounded Simplex-Structured Matrix Factorization", ICASSP 2022. [arXiv]

Deep orthogonal NMF

This Matlab code contains the hierarchical approach to tackle the deep ONMF problem, as well as other approaches and some data sets; see the paper 

P. De Handschutter and N. Gillis, "Deep orthogonal matrix factorization as a hierarchical clustering technique", EUSIPCO 2021, 23-27 August, Dublin. [pdf] 

BMME, block alternating Bregman Majorization Minimization with Extrapolation

This code  provides the implementation of BMME, a novel framework of block alternating Bregman majorization minimization for non-convex non-smooth optimization where the block function are not necessarily L-smooth, for two applications: orthogonal NMF, and the matrix completion problem; see the paper 

L. T. K. Hien, D.N. Phan, N. Gillis, M. Ahookhosh,  and P. Patrinos, "Block Alternating Bregman Majorization Minimization with Extrapolation", SIAM J. on Mathematics of Data Science 4 (1), pp. 1-25, 2022. [doi] [arXiv] 

Non-symmetric PSD Procrustes

This Matlab code contains several algorithms for solving the non-symmetric positive semidefinite (PSD) Procrustes problem: given rectangular matrices X and B, find the matrix A that minimizes the Frobenius norm of AX-B and such that A+AT is PSD. It contains our semi-analytical approach combining a clever initialization strategy and the fast gradient method. It also contains a function solving the problem with CVX. I can also solve the symmetric PSD problem (see below), and the complex non-PSD Procrustes problem; see 

M. K. Baghel, N. Gillis and P. Sharma, "On the non-symmetric semidefinite Procrustes problem", Linear Algebra and its Applications 648, pp. 133-159, 2022. [doi] [arXiv] 

Nonlinear subspace clustering

This code  contains several nonlinear subspace clustering algorithms, as well as data sets. It can be used to perform the numerical experiments of the following paper 

M. Abdolali and N. Gillis, "Beyond Linear Subspace Clustering: A Comparative Study of Nonlinear Manifold Clustering Algorithms", Computer Science Review 42, 100435, 2021. [doi] [arXiv] [Matlab] 

iADMM, an inertial alternating direction method of multipliers

This code  provides the implementation of iADMM, a novel framework of inertial alternating direction method of multipliers for non-convex non-smooth optimization, for an application in a low-rank decompositions problem; see the paper 

L. T. K. Hien,  D.N. Phan and N. Gillis, "Inertial alternating direction method of multipliers for non-convex non-smooth optimization", Computational Optimization and Applications 83, pp. 247-285, 2022. [doi] [arXiv]

Linear-quadratic NMF under separability

This code  provides the implementation of two new algorithms that allow to tackle linear-quadratic NMF that solve the problem X ≈ WH + (W ◦ W) H', where W,H,H' are nonnegative matrices and (W ◦ W) contains the component-wise products of any two columns of W. It relies on the assumption that each column of W appears purely in the data set (separability assumption); see the paper 

C. Kervazo, N. Gillis, and N. Dobigeon, "Provably robust blind source separation of linear-quadratic near-separable mixtures",  SIAM J. on Imaging Sciences 14 (4), pp. 1848-1889, 2021. [doi] [arXiv] 

Sparse Nonnegative Least Squares

This gitlab repository  provides the implementation of several methods to solve sparse nonnegative least squares problems; see the papers 

N. Nadisic, A. Vandaele, N. Gillis and J.E. Cohen, "Exact Sparse Nonnegative Least Squares",  ICASSP 2020, May 4-8, 2020, Barcelona. [pdf] [video]

N. Nadisic, A. Vandaele, J.E. Cohen and N. Gillis, "Sparse Separable Nonnegative Matrix Factorization", ECML-PKDD 2020. [arXiv] [pdf] [video]

N. Nadisic, J.E. Cohen, A. Vandaele, and N. Gillis, "Matrix-wise ℓ0-constrained Sparse Nonnegative Least Squares", Machine Learning, 2022. [doi] [arXiv] 

TITAN, an inertial block majorization-minimization algorithm

This code  provides the implementation of TITAN, a novel inerTial block majorIzation minimization framework for non-smooth non-convex opTimizAtioN problems, for two applications, namely sparse NMF and matrix completion; see the paper 

L. T. K. Hien,  D.N. Phan and N. Gillis, "An Inertial Block Majorization Minimization Framework for Nonsmooth Nonconvex Optimization", October 2020. [arXiv] 

Penalized β-NMF with disjoint linear constraints

This code  provides the implementation of multiplicative updates for three NMF models minimizing Dβ(V,WH) over variables W≥0 and H≥0, namely (1) β-NMF with sum-to-one constraints on the columns of H, (2) minimum-volume NMF with sum-to-one constraints on the columns of W, and (3) sparse NMF with ℓ1 penalty for the rows of H and ℓ2 norm csontraints on the columns of W. These three algorithms were designed following the framework described in the paper 

V. Leplat,  N. Gillis and J. Idier, "Multiplicative Updates for NMF with β-Divergences under Disjoint Equality Constraints", SIAM J. on Matrix Analysis and Applications 42 (2), 730-752, 2021. [doi] [arXiv]

Kullback-Leibler NMF

This code  provides various algorithms to tackle NMF using the Kullback-Leibler divergence as the objective function. This includes three new algorithms; see the paper 

L. T. K. Hien and N. Gillis, "Algorithms for Nonnegative Matrix Factorization with the Kullback-Leibler Divergence", Journal of Scientific Computing 87, 93, 2021. [doi] [arXiv] 

Deep MF

 This Matlab code allows you to run the two examples presented in the paper below, by tackling deep MF, decomposing a matrix X as W1W2...WLHL; see the paper 

P. De Handschutter, N. Gillis, and X. Siebert, "A survey on deep matrix factorizations", Computer Science Review 42, 100423, 2021. [doi] [arXiv]

Simplex-Structured Matrix Factorization (SSMF) via facet identification 

This Matlab code  allows you to tackle SSMF which is the following problem: Given an m-by-n matrix X and a factorization rank r, find an m-by-r matrix W and an r-by-n nonnegative matrix H with ||H(:,j)||1=1 for all j (that is, the columns of H belong to the unit simplex) such that X = WH. This code is based on the identification of the facets of conv(W), and works in the presence of noise and outliers; see the paper 

M. Abdolali and N. Gillis, "Simplex-Structured Matrix Factorization: Sparsity-based Identifiability and Provably Correct Algorithms",  SIAM J. on Mathematics of Data Science 3 (2), 593-623, 2021. [doi] [arXiv]

Minimax NMF

This Matlab code can be used to compute the  minimax NMF model: Given a set of matrices Xi (1≤i≤n) with the same number of rows, it comptes W≥0 and Hi≥0 to minimize max1≤i≤n ||Xi - WHi||F2 + λ logdet(WTW) where the second term is a volume regularizer; see the paper 

T. Marrinan and N. Gillis, "Hyperspectral Unmixing with Rare Endmembers via Minimax Nonnegative Matrix Factorization", EUSIPCO 2020. [pdf]  

Minimum enclosing ball of a collection of linear subspaces

This Matlab code can be used to compute the linear subspace whose maximal distance to a set of linear subspaces is minimized; see the paper 

T. Marrinan, P.-A. Absil and N. Gillis, "On a minimum enclosing ball of a collection of linear subspaces", Accepted in Linear Algebra and its Applications, 2021. [arXiv]

Off-diagonal Symmetric NMF

This Matlab code can be used to compute symmetric NMFs where the diagonal entries of the input matrix X are not taken into account. More precisely, it computes a nonnegative matrix H such that ∑i≠j |A - HHT|i,jp  is minimized where p = 1 or 2 (to the best of our knowledge, this is the first code that handles the case p=1); see the paper 

F. Moutier, A. Vandaele and N. Gillis, "Off-diagonal Symmetric Nonnegative Matrix Factorization", Numerical Algorithms, 2021. [doi] [arXiv] [pdf]

Grouped sparse projection

This Matlab code can be used to project a set of vectors {x1, x2,..., xn} of possible different dimensions onto a set of vectors {y1, y2,..., yn} such that (1) yi is close to xi for all i, and (2) the average sparsity if the yi's is at least s, where s ∈ [0,1] is a given parameter. The sparsity is defined using the measure of Hoyer based on the ratio between the ℓ1 and ℓ2 norms. This code can also be used to solve the NMF problem with sparsity constraints. A version in Python is also available, along with the numerical experiments on using this projection to design autoencoders;  see the paper 

N. Gillis, R. Ohib, S. Plis and V. Potluru, "Grouped Sparse Projection for Deep Learning", ICLR Worshop Hardware Aware Efficient Training (HAET), 2021. [arXiv] 

Near-Convex Archetypal Analysis (NCAA)

This Matlab code can be used to tackle a particular nonnegative matrix factorization (NMF) model closely related to archetype analysis (AA): Given a nonnegative data matrix X, a nonnegative matrix Y (e.g., Y=X) and an integer r, find nonnegative matrices A (with r columns) and H (with r rows) such that X≈YAH; see the paper

P. De Handschutter, N. Gillis, A. Vandaele and X. Siebert, "Near-Convex Archetypal Analysis", IEEE Signal Processing Letters 27 (1), pp. 81-85, 2020. [doi] [arXiv] [pdf] 

Robust Successive Projection Algorithm (RSPA)

This Matlab code can be used to tackle separable nonnegative matrix factorization (separable NMF): Given a nonnegative data matrix X and an integer r, find an index set K with r indices such that there exists a nonnegative matrix H with X≈X(:,K)H; see the paper 

N. Gillis, "Successive Projection Algorithm Robust to Outliers", CAMSAP 2019, December 15-18 2019, Guadeloupe. [arXiv]

Minimal-norm static feedbacks using dissipative Hamiltonian matrices

This Matlab code can be used to compute minimal-norm static-state and static-output feedback matrices. More precisely, given a continuous LTI system of the form x'(t) = Ax(t) + Bu(t), y(t) = Cx(t), this code tackles the problem of finding a matrix K with minimal norm such that A-BK is stable (static-state feedback problem) or such that A-BKC is stable (static-output feedback problem); see the paper 

N. Gillis, P. Sharma, "Minimal-norm static feedbacks using dissipative Hamiltonian matrices", Linear Algebra and its Applications 623, pp. 258-281, Special issue in honor of Paul Van Dooren, 2021. [doi] [arXiv] [pdf]

Minimum-Volume β-NMF

This Matlab code tackles the following optimization problem: Given an m-by-n nonnegative matrix V and a factorization rank r, min-vol β-NMF is the problem of finding an m-by-r nonnegative matrix W and an r-by-n nonnegative matrix H with ||W(:,j)||1=1 for all j such that Dβ(V,WH) + λ logdet(WTW+ δ I) is minimized, where Dβ(V,WH) is the β-divergence between V and WH. The code handles β=1 (Kullback-Leibler divergence) and β=0 (Itakura-Saito divergence); see the paper 

V. Leplat, N. Gillis, A.M.S. Ang, "Blind Audio Source Separation with Minimum-Volume Beta-Divergence NMF", IEEE Transactions on Signal Processing 68, pp. 3400-3410, 2020. [doi] [arXiv] 

Generalized separable NMF

This Matlab code tackles the generalized separable nonnegative matrix factorization (GS-NMF) problem: Given a nonnegative matrix M and an integer r, find index sets K1 and K2 with a total of r indices, and nonnegative matrices P1 and P2, such that || M - M(:,K1) P1 + P2 M(K2,:) ||F is minimized; see the paper 

J. Pan and N. Gillis, "Generalized Separable Nonnegative Matrix Factorization", IEEE Trans. on Pattern Analysis and Machine Intelligence 43 (5), pp. 1546-1561, 2021. [doi] [arXiv] [pdf]

Inertial Block Mirror Descent Methods for NMF

This Matlab code allows you to solve the following nonnegative matrix factorization (NMF) problem: Given X and r, compute W≥0 and H ≥0 such that ||X-WH||F is minimized. To achieve this goal, the code uses the inertial block mirror descent methods described in the paper 

L. T. K. Hien, N. Gillis, P. Patrinos, "Inertial Block Proximal Methods for Non-Convex Non-Smooth Optimization", ICML 2020. [arXiv] 

Distributionally Robust and Multi-Objective NMF

This Matlab code allows you to solve the NMF problem where the objective function is a weighted sum of several β-divergences measuring the error between the input matrix X and the factorization WH, where W≥0 and H≥0. It also allows you to solve the distributionally robust NMF problem where the goal is to minimize the maximum value among the β-divergences; see 

N. Gillis, L. T. K. Hien, V. Leplat and V. Y. F. Tan, "Distributionally Robust and Multi-Objective Nonnegative Matrix Factorization", IEEE Trans. on Pattern Analysis and Machine Intelligence, 2021. [doi] [arXiv]

Approximating the nearest Ω-Stable Matrix

This Matlab code can be used to compute the nearest Ω-stable matrix from an unstable one, that is, to compute a matrix whose eigenvalues are within Ω and that is as close as possible to the given input matrix. Our algorithm uses a block coordinate descent method applied on a reformulation of the problem. The set Ω, within the complex plane, is the intersections of conic sectors, vertical strips and disks. 

See N. Choudhary, N. Gillis, P. Sharma, "On approximating the nearest Ω-stable matrix", Numerical Linear Algebra with Applications 27 (3), e2282, 2020. [doi] [arXiv] [pdf] 

Minimum-Volume NMF

This Matlab code tackles the following optimization problem: Given an m-by-n nonnegative matrix X and a factorization rank r, min-vol NMF is the problem of finding an m-by-r nonnegative matrix W and an r-by-n nonnegative matrix H with ||H(:,j)||1≤1 for all j such that ||X - WH||2F + λ logdet(WTW+ δ I) is minimized. 

See V. Leplat, A.M.S. Ang, N. Gillis, "Minimum-Volume Rank-Deficient Nonnegative Matrix Factorizations", ICASSP 2019, May 12-17, 2019, Brighton, UK. [pdf] 

Identifiability of Complete Dictionary Learning

This Matlab code can be used to generate two low-rank sparse SCA decompositions of the matrix M of dimensions r x (r3-2r2). The decompositions have the form M  = DiBi where rank(Di)=r, ||Bi(:,j)||0 = r-1 for all j, and each subspace spanned by r-1 columns of Di contains r2-2r columns of M whose spark is r, for i = 1,2. 

J.E. Cohen and N. Gillis, "Identifiability of Complete Dictionary Learning", SIAM J. on Mathematics of Data Science 1 (3), pp. 518-536, 2019. [doi] [arXiv] [pdf]

Approximating the nearest stable matrix pair in the discrete-time case

This Matlab code can be used to compute the nearest stable matrix pair from an unstable one in the discrete-time case. Our algorithm uses a block coordinate descent method applied on a reformulation of the problem. 

See N. Gillis, M. Karow and P. Sharma, "A note on approximating the nearest stable discrete-time descriptor system with fixed rank", Applied Numerical Mathematics 148, pp. 131-139, 2020. [doi] [arXiv] [pdf]

Nonnegative SVD with low-rank correction (NNSVD-LRC)

This matlab code contains a new SVD-based initialization for NMF, namely NNSVD-LRC. Compared to previous SVD-based initializations, it is faster and reduces the error as the factorization rank increases. 

See Atif M. Syed, Sameer Qazi, Nicolas Gillis, "Improved SVD-based Initialization for Nonnegative Matrix Factorization using Low-Rank Correction", Pattern Recognition Letters 122, pp. 53-59, 2019. [doi] [arXiv] [pdf]

Accelerating NMF algorithms using extrapolation

Given a nonnegative m-by-n matrix X and a factorization rank r, nonnegative matrix factorization (NMF) looks for two nonnegative matrices W and H of dimension m-by-r and r-by-n respectively such that X ≈ WH. This matlab code contains the acceleration of NMF algorithms described in the paper below. The idea is to use extrapolation, that is, use the update W* = Wn + β (Wn - W) where Wn and W are the previous two iterates and β ∈ [0,1] is a parameter (and similarly for H). 

See A.M.S. Ang and N. Gillis, "Accelerating Nonnegative Matrix Factorization Algorithms using Extrapolation", Neural Computation 31 (2), pp. 417-439, 2019. [doi] [arXiv] [pdf] [Matlab]

Stabilizing discrete-time linear systems

This Matlab code can be used to compute the nearest stable matrix from an unstable one in the discrete-time case. Our algorithm uses a fast gradient method applied on a reformulation of the problem: the square matrix A is stable if and only if there exists a symmetric positive definite matrix S, an orthogonal matrix U and a symmetric positive semidefinite matrix B with ||B||≤1 such that A = S-1 U B S. 

See N. Gillis, M. Karow and P. Sharma, "Stabilizing discrete-time linear systems", Linear Algebra and its Applications 573, pp. 37-53, 2019. [doi] [arXiv] [pdf] 

Scalable and Robust Sparse Subspace Clustering (SR-SSC)

This Matlab code can be used to perform the experiments described in the paper below. In particular, it can be used to run SR-SSC that performs subspace clustering using a randomized hierarchical clustering strategy and mutlilayer graphs. 

See M. Abdolali, N. Gillis and M. Rahmati, "Scalable and Robust Sparse Subspace Clustering Using Randomized Clustering and Multilayer Graphs", Signal Processing 163, pp. 166-180, 2019. [doi] [arXiv] [pdf]

Nearest Positive-Real System

This Matlab code contains a method to solve the following problem: given a linear-time invariant dynamical system defined by Ex˙ = Ax+Bu and y = Cx+Du, find a nearby positive-real system. Our method is based on a reformulation of this problem using port-Hamiltonian (PH) systems for which A=(J-R)Q, B=F-P, C=(F+P)T and D=S+N, where QTE and [R P; PT S] are positive semidefinite, N=-NT, and J=-JT. 

See N. Gillis and P. Sharma, "Finding the nearest positive-real system", SIAM J. on Numerical Analysis 56 (2), 1022-1047, 2018. [doi] [arXiv]

Low-Rank Matrix Approximation in the Infinity Norm

This Matlab code implements an exact cyclic coordinate descent method for the component-wise ℓ∞-norm low-rank matrix approximation problem: Given an m-by-n matrix M and a factorization rank r, find an m-by-r matrix U and an r-by-n matrix V such that ||M-UV||∞ = maxi,j |M-UV|ij is minimized. By default, it initializes the algorithm with the optimal solution of the ℓ2-norm problem using the truncated singular value decomposition provided by Matlab. 

See N. Gillis and Y. Shitov, "Low-Rank Matrix Approximation in the Infinity Norm", Linear Algebra and its Applications 581, pp. 367-382, 2019. [doi] [arXiv] [pdf] 

Nearest Stable Matrix Pairs

This Matlab code contains the fast gradient method to solve the following problem: given an unstable matrix pair (E,A), find the stable matrix pair (M,X) such that ||A-X||²F + ||E-M||²F is minimized. The algorithm is based on the following reformulation: X = (J-R)Q where J is antisymmetric (J=-JT), R is positive semidefinite, Q is invertible and QTM is positive semidefinite. The pair ((J-R)Q,M) is a so-called dissipative Hamiltonian (DH) matrix pair.  

See N. Gillis, V. Mehrmann and P. Sharma, "Computing nearest stable matrix pairs", Numerical Linear Algebra with Applications 25 (5), e2153, 2018. [doi] [arXiv] [pdf]

PSD Procrustes

This Matlab code contains several algorithms for solving the positive semidefinite Procrustes problem: given rectangular matrices X and B, find the symmetric positive semidefinite matrix A that minimizes the Frobenius norm of AX-B. It contains the projected gradient method, the fast gradient method, a method from the paper 'Andersson & Elfving, A constrained Procrustes problem, SIAM Journal on Matrix Analysis and Applications, 18(1), 124-139, 1997', and our semi-analytical approach combining a clever initialization strategy and the fast gradient method. It also contains a function solving the problem with CVX. 

See N. Gillis and P. Sharma, "A semi-analytical approach for the positive semidefinite Procrustes problem", Linear Algebra and its Applications 540, pp. 112-137, 2018. [doi] [arXiv]

Nearest Stable Matrix

 This Matlab code provides the codes for three algorithms to solve the following problem: given an n-by-n matrix A, find the n-by-n matrix X such that ||A-X||F is minimized and X is stable. These algoritms are based on the following reformulation: X = (J-R)Q where J is antisymmetric (J=-JT), and R and Q a positive semidefinite. 

This code also contains three algorithms kindly provided to us by their authors: 

- SuccConv from the paper 'F.-X. Orbandexivry, Yu. Nesterov, and P. Van Dooren, Nearest stable system using successive convex approximations, Automatica, 49 (2013), pp. 1195-1203',  

- BFGS from Michael Overton using his toolbox HANSO, and

- A method based on GRANSO that solves non-smooth optimization problems; see the paper here. The use of GRANSO for the nearest stable matrix problem was implemented by Tim Mitchell. Note that comparisons with this algorithm are not provided in our paper since Tim gave us this algorithm later on. For large matrices, our code still works best (for small matrices, it depends to which local minima each algorithm converges to; cf. the discussion in the paper). 

If you use these codes, please refer to the work of these authors. 

See N. Gillis and P. Sharma, "On computing the distance to stability for matrices using linear dissipative Hamiltonian systems", Automatica 85, pp. 113-121, 2017. [doi] [arXiv]

Nonnegative Sparse Regression with Self Dictionary

This Matlab code implements a fast gradient method to solve the following problem: given an m-by-n matrix M, find an n-by-n nonnegative matrix X such that M ≈ MX and such that X has a few non-zero rows. The matrix X allows to recover a subset of the columns of M (corresponding to the non-zero rows of X) whose convex cone contains approximately all the columns of M (since X is nonnegative). In hyperspectral imaging, this method allows to identify pure pixels. 

See N. Gillis and R. Luce, "A Fast Gradient Method for Nonnegative Sparse Regression with Self Dictionary", IEEE Trans. on Image Processing 27 (1), pp. 24-37, 2018. [doi] [arXiv]

Component-Wise ℓ1-Norm Low-Rank Matrix Approximation

This Matlab code implements an exact cyclic coordinate descent method for the component-wise ℓ1-norm low- rank matrix approximation problem: Given an m-by-n matrix M and a factorization rank r, find an m-by-r matrix U and an r-by-n matrix V such that ||M-UV||1 = sumi,j |M-UV|ij is minimized. By default, it initializes the algorithm with the optimal solution of the ℓ2-norm problem using the truncated singular value decomposition provided by Matlab. 

See N. Gillis and S.A. Vavasis, On the Complexity of Robust PCA and ℓ1-Norm Low-Rank Matrix Approximation, Mathematics of Operations Research 43 (4), pp. 1072-1084, 2018. [arXiv] [pdf] [doi]

Symmetric NMF

This Matlab code solves symmetric nonnegative matrix factorization (symNMF) problems using an exact coordinate descent method. Given an n-by-n symmetric nonnegative matrix A and a factorization rank r, symNMF is the problem of finding an n-by-r nonnegative matrix H such that ||A - HHT||F is minimized.  

See A. Vandaele, N. Gillis, Q. Lei, K. Zhong, I. Dhillon, Coordinate Descent Methods for Symmetric Nonnegative Matrix Factorization, IEEE Transactions on Signal Processing 64 (21), pp. 5571-5584, 2016. [doi] [arXiv]

Clique Finding Algorithm

This Matlab code tries to identify a maximum clique in a graph (the input is the adjacency matrix of that graph). The preprocessed la1 data set can be downloaded here. 

See M.T. Belachew and N. Gillis, "Solving the Maximum Clique Problem with Symmetric Rank-One Nonnegative Matrix Approximation", Journal of Optimization Theory and Applications, 173 (1), pp. 279-296, 2017. [doi] [arXiv]

Heuristics for Exact NMF

The code available on this website tries to compute exact NMF's: given an m-by-n nonnegative matrix M and a factorization rank r, it looks for an m-by-r nonnegative matrix U and an r-by-n nonnegative matrix V such that M = UV; see the website for all the details, and the paper A. Vandaele, N. Gillis, F. Glineur and D. Tuyttens, "Heuristics for Exact Nonnegative Matrix Factorization", Journal of Global Optimization 65 (2), pp 369-400, 2016. [doi] [arXiv] 

Exact and Heuristic Algorithms for Semi-NMF

Given an m-by-n matrix M, semi-NMF looks for an m-by-r matrix U and a r-by-n nonnegative matrix V such that M ≈ UV. This matlab code contains: 

- A method that computes the semi-nonnegative rank (smallest r such that M = UV under the conditions above) and a corresponding factorization M = UV with V≥0. 

- The initialization strategy based on k-means proposed in 'Convex and Semi-Nonnegative Matrix Factorizations', Ding, Li and Jordan, IEEE Trans. on Pattern Analysis and Machine Intelligence 32 (1), 2010. 

- Two initialization strategies for semi-NMF described in the paper below. 

- File to run experiments on synthetic data sets to compare the different initialization strategies. 

See N. Gillis and A. Kumar, "Exact and Heuristic Algorithms for Semi-Nonnegative Matrix Factorization", SIAM J. on Matrix Analysis and Applications 36 (4), pp. 1404-1424, 2015. [doi] [arXiv] [Slides]

Enhancing Pure-Pixel Identification via Preconditioning

This matlab code contains: 

- Several pure-pixel search algorithms: VCA, SPA, XRAY (a.k.a. near-separable NMF algorithms).  

- Experiments on synthetic data sets comparing the different preconditionings (namely, the SDP-based preconditioning, pre-whitening, and the SPA-based preconditioning). 

See N. Gillis and W.-K. Ma, "Enhancing Pure-Pixel Identification Performance via Preconditioning", SIAM J. on Imaging Sciences 8 (2), pp. 1161-1186, 2015. [arXiv] [Slides] [doi] [pdf]

Successive Nonnegative Projection Algorithm (SNPA)

This matlab code contains: 

- The successive projection algorithm (SPA); see arXiv:1208.1237, and below.

- The successive nonnegative projection algorithm (SNPA). 

- An implementation of the fast conical hull algorithm of Kumar, Sindhwani & Kambadur (ICML '13); see arXiv:1210.1190 (XRAY).

- Examples on synthetic data sets (Rank-deficient/Ill-conditioned and Dirichlet/Middle points experiments; cf. v3 of arXiv paper).

See N. Gillis, "Successive Nonnegative Projection Algorithm for Robust Nonnegative Blind Source Separation", SIAM J. on Imaging Sciences 7 (2), pp. 1420-1450, 2014. [doi] [arXiv] [pdf]

Remark.  A better implementation of SNPA that runs in O(nnz * r) operations, where nnz is the number of non-zeros of the input matrix and r is the factorization rank, is available here; see the discussion the NMF book (pages 235-236). 

SDP Based Preconditioning for More Robust Near-Separable NMF

This matlab code contains: 

- The successive projection algorithm (SPA); see arXiv:1208.1237, and below.

- The post-processed SPA; see arXiv:1212.4777 (Post-SPA).

- The preconditioned SPA (Prec-SPA). 

- The post-processed preconditioned SPA (Post-Prec-SPA).

- The vertex component analysis (VCA) of Nascimento and Bioucas-Dias.

- An implementation of the fast conical hull algorithm of Kumar, Sindhwani & Kambadur (ICML '13); see arXiv:1210.1190 (XRAY).

- Two examples on synthetic data sets (Middle points and Hubble; see paper below).

You can find here the Hubble telescope hyperspectral image (kindly provided to us by Robert J. Plemmons). 

See N. Gillis and  S.A. Vavasis, "Semidefinite Programming Based Preconditioning for More Robust Near-Separable Nonnegative Matrix Factorization", SIAM J. on Optimization 25 (1), pp. 677-698, 2015. [arXiv] [Slides] [doi] [pdf]

Hierarchical Clustering of Hyperspectral Images with NMF

This matlab code provides a hierarchical clustering algorithm particularly efficient for high-resolution hyperspectral images. It is based on rank-two NMF and convex geometry. Based on the clustering, it also provides the spectral signatures of pure pixels in such images. 

You can find here the Urban hyperspectral image (as a 307x307x162 tensor). 

See N. Gillis, D. Kuang and H. Park, "Hierarchical Clustering of Hyperspectral Images using Rank-Two Nonnegative Matrix Factorization", IEEE Trans. on Geoscience and Remote Sensing 53 (4), pp. 2066-2078, 2015. [arXiv] [doi]

Robust near-separable NMF using LP

This matlab code contains: 

- A CPLEX code for HottTopixx of Bittorf, Recht, Ré, and Tropp (NIPS '12); see arXiv:1206.1270. (See below for a CVX version.)

- An implementation of the fast and recursive robust algorithm of Gillis and Vavasis; see arXiv:1208.1237, and below.

- An implementation of the fast conical hull algorithm of Kumar, Sindhwani & Kambadur (ICML '13); see arXiv:1210.1190.

- A CPLEX code for our new LP model for near-separable NMF; see paper below.

- A test file minibench.m to compare the above algorithms on synthetic datasets.

- A test file Swimmer_test.m to compare the above algorithm on the swimmer data set. 

See N. Gillis and R. Luce, "Robust Near-Separable Nonnegative Matrix Factorization Using Linear Optimization", Journal of Machine Learning Research 15 (Apr), pp. 1249-1280, 2014. [pdf] [Slides] [Video]

Robustness analysis of Hottopixx, an LP model for NMF

This matlab code provides the constructions decribed in the paper below. It also contains a CVX code for the LP model from Bittorf, Recht, Ré, and Tropp, "Factoring nonnegative matrices with linear programs", Advances in Neural Information Processing Systems 25, pp. 1223-1231, 2012. (Note: Third version of the code available, there was an error in the post-processing.) 

See N. Gillis, "Robustness Analysis of Hottopixx, a Linear Programming Model for Factoring Nonnegative Matrices", SIAM J. Matrix Anal. & Appl. 34 (3), pp. 1189-1212, 2013. [doi] [arXiv]

Fast and robust recursive algorithm for separable NMF

A nonnegative matrix M = UV with U≥0 and V≥0 is separable if each column of U is equal to a column of M up to a scaling factor. In other words, there exists a cone spanned by a small subset of the columns of the input nonnegative data matrix M containing all columns. This matlab code provides an implementation for extracting, among the columns of M, the ones corresponding to the columns of U. This algorithm is a generalization of the successive projection algorithm (SPA). 

See N. Gillis and  S.A. Vavasis, "Fast and Robust Recursive Algorithms for Separable Nonnegative Matrix Factorization", IEEE Trans. on Pattern Analysis and Machine Intelligence 36 (4), pp. 698-714, 2014.  [doi] [arXiv] [Slides]

Biclique finding algorithm

Given the biadjacency matrix of a bipartite graph, this matlab code returns a large biclique of that graph (that is, a complete bipartite induced subgraph). It also contains two running examples: one with synthetic data, and one with the text dataset la1.

See N. Gillis and F. Glineur, "A Continuous Characterization of the Maximum-Edge Biclique Problem”, Journal of Global Optimization 58 (3), pp. 439-464, 2014. [doi] [pdf]

Accelerated MU and HALS algorithms for NMF

Given a nonnegative m-by-n matrix M, NMF looks for two nonnegative matrices U and V of dimension m-by-r and r-by-n respectively such that M ≈ UV. This matlab code provides an implementation for the following methods:

- Accelerated multiplicative updates of Lee and Seung (A-MU).

- Accelerated hierarchical alternating least squares of Cichocki et al. (A-HALS).

- Accelerated projected gradient method of Lin (A-PG).

- A comparison of the different algorithms on the popular CBCL face database (RunMe.m).

All methods aim at minimizing the Frobenius norm (that is, sum of the squared entries) of the residual ||M-UV||F. A-HALS is the most efficient method (on the tested datasets). 

See N. Gillis and F. Glineur, "Accelerated Multiplicative Updates and Hierarchical ALS Algorithms for Nonnegative Matrix Factorization”, Neural Computation 24 (4), pp. 1085-1105, 2012. [doi] [arXiv]

Preprocessing for sparse and unique NMF

Given a nonnegative matrix M (where no column of M is multiple of another column), this matlab code computes an inverse positive matrix Q (that is, the inverse of Q exists and is nonnegative) so that P(M) = MQ is nonnegative and sparse. This preprocessing allows to obtain a more well-posed NMF problem whose solutions are sparser. We have that any NMF (U,V') of P(M) provides an NMF for M since M=P(M)Q-1≈U(V'Q-1)=UV where V=(V'Q-1)≥0.  

This code also contains a very efficient sparse NMF method based on the accelerated HALS algorithm (see above) and l1-norm penalty. 

See N. Gillis, "Sparse and Unique Nonnegative Matrix Factorization Through Data Preprocessing", Journal of Machine Learning Research 13 (Nov), pp. 3349-3386, 2012. [pdf] [Slides]

Geometric interpretation of the nonnegative rank

This matlab code contains two main methods:

- Given a nonnegative rank-three matrix M, NMFrank3.m displays the 2-D geometric interpretation of the exact NMF problem for M. 

- Given a nonnegative rank-three matrix M, and a rank four nonnegative matrix U so that there exists a nonnegative matrix V with M=UV, Visualisation3D.m displays the 3-D geometric interpretation of the corresponding exact NMF problem. 

See N. Gillis and F. Glineur, "On the Geometric Interpretation of the Nonnegative Rank ”, Linear Algebra and its Applications 437 (11), pp. 2685-2712, 2012. [doi] [arXiv] [Slides]

Nonnegative matrix underapproximation (NMU)

a) Global NMU solves NMF problems with the additional underapproximation constraint UV ≤ M which allows to obtain better part-based decompositions. 

See N. Gillis and F. Glineur, "Using Underapproximations for Sparse Nonnegative Matrix Factorization", Pattern Recognition 43 (4), pp. 1676-1687, 2010. [doi] [arXiv]

b) Recursive NMU solves the NMF problem recursively using rank-one NMU, computing one rank-one factor at a time. 

See paper above, and N. Gillis and R.J. Plemmons, "Dimensionality Reduction, Classification, and Spectral Mixture Analysis using Nonnegative Underapproximation", Optical Engineering 50, 027001, February 2011. [doi] [pdf] [Slides]

c) Sparse NMU implements recursive NMU with additional sparsity constraints on the factor U. It makes it more efficient to extracting sparse abundance maps in hyperspectral images.  

See N. Gillis and R.J. Plemmons, "Sparse Nonnegative Matrix Underapproximation and its Application for Hyperspectral Image Analysis", Linear Algebra and its Applications 438 (10), pp. 3991-4007, 2013. [doi] [Slides]

d) Spatial NMU implements recursive NMU with additional spatial constraints on the factor U, that is, it takes into account the fact that neighboring pixels are more likely to contain the same materials. 

See N. Gillis, R.J. Plemmons and Q. Zhang, "Priors in Sparse Recursive Decompositions of Hyperspectral Images", in Algorithms and Technologies for Multispectral, Hyperspectral, and Ultraspectral Imagery XVIII, edited by Sylvia S. Shen, Paul E. Lewis, Proceedings of SPIE Vol. 8390 (2012). [doi] [pdf]

e) Prior NMU implements recursive NMU with additional spatial and sparsity constraints on the factor U, that is, it takes into account the fact that neighboring pixels are more likely to contain the same materials and the fact that features should contain a relatively small number of pixels. 

See G. Casalino and N. Gillis, "Sequential Dimensionality Reduction for Extracting Localized Features", Pattern Recognition 63, pp. 15-29, 2017. [doi] [arXiv]