nonnegative matrix factorization, polyhedral cones, double description method, obtuseness, alternating linear programming, exact factorization

1 Introduction↩︎

Nonnegative matrix factorisation seeks \(W \in \mathbb{R}^{m \times r}_{\ge 0}\) and \(H \in \mathbb{R}^{n \times r}_{\ge 0}\) such that \(A \approx W H^{\!\top}\) for a given \(A \in \mathbb{R}^{m \times n}_{\ge 0}\) of nonnegative rank at most \(r\). Two distinct problems hide inside this template. The approximate problem, used in topic modelling, image decomposition and recommendation systems, is typically attacked by multiplicative updates [1], hierarchical alternating least squares [2], or block coordinate descent [3], and is NP-hard in the worst case [4]. The exact problem assumes the data \(A\) admits an \(r\)-factorisation with zero residual and asks for a certified \((W, H)\) recovering it; structural conditions such as separability [5]–[8] make this tractable. This paper is firmly on the exact side, but without imposing separability.

We develop a small, fully executable toolkit (svd.wl, svdGb.wl, svd_obtuse.wl) in Wolfram Language that combines four classical ingredients into one pipeline:

1. a truncated SVD \(A = U_r S_r V_r^{\!\top}\),

2. the polyhedral cone \(\mathop{\mathrm{cone}}(A_{\mathrm{orth}})\) of nonnegative preimages and its dual extreme-ray representation via the Double Description Method [9], [10] and cddlib,

3. a uniform-support parameterisation of \(Q, P\) on \(r\)-subsets \(T, K\) of cone-ray indices whose joint constraint reduces to the entrywise nonnegativity of an \(r \times r\) witness matrix \(M_T\) (Proposition [prop:Mt]), and

4. when \(M_T\) is infeasible, an augmented \(r \times (r + k)\) slack-LP alternation that recovers a witness on enlarged supports.

The output is a certificate, not an estimate: when the toolkit succeeds it returns \((W, H)\) such that \(\|A - W H^{\!\top}\|_F / \|A\|_F\) is at machine precision.

1.0.0.1 Contributions.

The technical contributions of this paper are:

1. Closed-form witness under uniform support. For exact-rank inputs, the joint nonnegativity constraint \(Q P^{\!\top} = I_r\) collapses to the entrywise condition \(M_T = R_T^{-1} (R_K^{\!\top})^{-1} \ge 0\) for an \(r\)-subset pair \((T, K)\). This is a microsecond test once \(R\) is known.

2. Obtuseness-ranked support search and its saturation. Ranking \(r\)-subsets by their geometric near-orthogonality and walking the top of each list extends reach from \(r = 4\) (lex enumeration) to \(r = 4, 5, 6\) at \(m = n = 10\). We characterise the saturation curve in \(100\)-trial Monte Carlo at three budgets (top-\(5\), top-\(200\), top-\(400\)), and show the heuristic saturates at top-\(200\) and that the residual failures are structural rather than budget-limited.

3. Scaling walls. We document where the pipeline breaks: at fixed \(r\) the success rate decreases as \(m, n\) grows (a \(100\)-trial sweep at \(m = n = 15\) collapses to \(37/7/0/0/0\) at \(r = 4, \dots, 8\)), and on Olivetti (\(400 \times 4096\)) cddlib itself does not finish in a one-hour budget on the \(n\)-side even at \(r = 8\).

4. Augmented-support alt-LP fallback (new). At each pair where the closed-form \(M_T\) is infeasible we augment both supports by \(k = 2\) rays maximally angularly-separated from the existing ones and solve a short alt-LP on the resulting \(r \times (r + k)\) bases. This raises the same \(m = n = 10\) success counts from \(79/85/58\) to \(\mathbf{99/95/75}\) at \(r = 4, 5, 6\) — breaking the structural ceiling that top-\(400\) alone cannot.

The executable toolkit and JSON results of every Monte Carlo are released alongside the paper.

2 Preliminaries↩︎

2.0.0.1 SVD truncation.

For \(A \in \mathbb{R}^{m \times n}_{\ge 0}\) of nonnegative rank \(r\) we compute the thin SVD \(A = U S V^{\!\top}\) and form \[A_{\mathrm{orth}}= U_{:,1:r} \sqrt{S_{1:r, 1:r}}, \qquad A_{\mathrm{orth}1}= V_{:,1:r} \sqrt{S_{1:r, 1:r}}, \label{eq:Aorth}\tag{1}\] so that \(A = A_{\mathrm{orth}}A_{\mathrm{orth}1}^{\!\top}\) and any factorisation \(W = A_{\mathrm{orth}}Q\), \(H = A_{\mathrm{orth}1}P\) with \(Q P^{\!\top} = I_r\) reconstructs \(A\) exactly.

2.0.0.2 Cones and the double description method.

The set of \(x \in \mathbb{R}^r\) with \(A_{\mathrm{orth}}x \ge 0\) is a polyhedral cone \(\mathop{\mathrm{cone}}(A_{\mathrm{orth}}) \subset \mathbb{R}^r\). Its extreme rays form an \(r \times k_1\) matrix \(R_{\mathrm{orth}}\) obtainable from the \(m \times r\) matrix of facet normals by the Double Description Method (DDM) [9], [10]. Define \(R_{\mathrm{orth}1}\in \mathbb{R}^{r \times k_2}\) analogously for \(A_{\mathrm{orth}1}\). Any column of \(Q\) (resp. \(P\)) that yields nonnegative \(W\) (resp. \(H\)) must lie in \(\mathop{\mathrm{cone}}(R_{\mathrm{orth}})\) (resp. \(\mathop{\mathrm{cone}}(R_{\mathrm{orth}1})\)). We rely on cddlib [10] for the DDM and on Mathematica’s LinearProgramming, GroebnerBasis and Reduce [11], [12] for the LP and algebraic sub-problems.

2.0.0.3 Nonneg LPs and slack relaxation.

Throughout we use the standard slack-LP relaxation of a feasibility problem \(A x = b, \; x \ge 0\): introduce \(s_+, s_- \ge 0\) and minimise \(\mathbf{1}^{\!\top}(s_+ + s_-)\) subject to \(A x + s_+ - s_- = b\) and \(x, s_\pm \ge 0\). The optimum is zero iff the original system is feasible.

3 The Uniform-Support Witness↩︎

Pick an \(r\)-subset \(T \subset [k_1]\) of cone-ray indices on the \(A_{\mathrm{orth}}\)-side and an \(r\)-subset \(K \subset [k_2]\) on the \(A_{\mathrm{orth}1}\)-side, and restrict \(Q, P\) to be supported on \(T, K\): \(Q = R_T \mu\) and \(P = R_K \nu\) with \(\mu, \nu \in \mathbb{R}^{r \times r}_{\ge 0}\). The joint condition \(Q P^{\!\top} = I_r\) becomes \[R_T \, \mu \, \nu^{\!\top} R_K^{\!\top} = I_r. \label{eq:joint}\tag{2}\]

Proposition [prop:Mt] is the engine of the toolkit. Given \(R_{\mathrm{orth}}, R_{\mathrm{orth}1}\) from the DDM step and a candidate index pair \((T, K)\), it takes one \(r \times r\) matrix inverse and one nonnegativity check to decide feasibility; if feasible the factorisation \((W, H) = (A_{\mathrm{orth}}R_T M_T, \; A_{\mathrm{orth}1}R_K)\) reconstructs \(A\) at machine precision (Section 5).

4 Obtuseness-Ranked Support Search↩︎

Because the candidate space \(\binom{k_1}{r} \times \binom{k_2}{r}\) quickly outgrows any enumeration budget, we must rank candidates. The geometric intuition behind Proposition [prop:Mt] is that \(M_T \ge 0\) is more likely for fat (near-orthogonal) sub-cones of \(R_{\mathrm{orth}}, R_{\mathrm{orth}1}\). We capture this with the obtuseness score: \[\mathrm{obtuseness}(T) \;:=\; \frac{|\det(R_T)|}{\prod_{i \in T} \|R_T[:, i]\|} \;\in\; [0, 1], \label{eq:obtuseness}\tag{4}\] which is \(1\) iff \(T\) is orthonormal and \(0\) iff \(T\) is rank-deficient.

The default \(\mathrm{maxTries} = 200\) admits at most \(40{,}000\) candidate pairs per instance; we will sweep this budget below.

5 Monte Carlo Characterisation at \(m = n = 10\)↩︎

To turn the qualitative “the heuristic works at small \(r\)” into distributional evidence, we run a \(100\)-trial Monte Carlo at \(m = n = 10\), \(r \in \{4, 5, 6\}\). Each trial draws \(W, H \sim \mathrm{U}[0, 1]^{10 \times r}\) and sets \(A = W H^{\!\top}\) (exact rank-\(r\), nonnegative). Per-trial seeds \(1000 r + \mathrm{trial}\) ensure reproducibility.

5.1 Baseline at \(\mathrm{maxTries} = 200\)↩︎

Table 1 reports the Monte Carlo result with the default budget. Two baselines (multiplicative updates [1] and HALS [2], both \(500\) iterations) are run on the same per-trial seeds for comparison.

5.2 Saturation curve: top-\(5\), top-\(200\), top-\(400\)↩︎

Whether the obtuseness heuristic is search-budget limited or structurally limited is settled by sweeping \(\mathrm{maxTries}\). Table 2 reports the three-point sweep.

5.2.0.1 Reading.

The saturation message is unambiguous: top-\(200\) is at the structural ceiling of what the closed-form \(M_T\) test can recover at this size. The remaining failures are not pairs hiding deep in the ranking but pairs whose cones intrinsically lack any near-orthogonal \(r\)-subset on at least one side (median best obtuseness in no-feasible trials drops monotonically from \(0.74\) at \(r=4\) to \(0.33\) at \(r=6\)). Pushing the budget further is not the right lever; lifting the uniform-support restriction is. We come back to this in Section 7.

6 Where the Pipeline Breaks↩︎

6.1 Larger \(m, n\) at fixed \(r\)↩︎

A natural conjecture is that \(m, n \gg r\) should help the cone toolkit: more rows in \(A\) give more facets of the cone and hence more rays from which to draw a near-orthogonal sub-basis. We test this with a parallel \(100\)-trial Monte Carlo at \(m = n = 15\), \(r \in \{4, \dots, 8\}\) (otherwise identical parameters to Table 1).

6.2 Olivetti faces (\(400 \times 4096\)): DDM stops working↩︎

We tried the same pipeline on the canonical NMF benchmark, the Olivetti faces dataset (\(A \in \mathbb{R}^{400 \times 4096}_{\ge 0}\)) at \(r = 25\) and \(r = 8\). The toolkit never reaches the obtuseness step — the prerequisite DDM on the \(n = 4096\) side exceeds budget even at the smallest \(r\) we tried:

• At \(r = 25\), \(\mathrm{ddmBudget} = 600\,\mathrm{s}\): cdd on \(A_{\mathrm{orth}}\) (a \(400 \times 25\) inequality system) times out; cdd on \(A_{\mathrm{orth}1}\) (\(4096 \times 25\)) is even larger.

• At \(r = 8\), \(\mathrm{ddmBudget} = 3{,}600\,\mathrm{s}\) (\(1\) hour) per cdd call: cdd on \(A_{\mathrm{orth}}\) (\(400 \times 8\)) finishes in \(\sim\!3\) minutes; cdd on \(A_{\mathrm{orth}1}\) (\(4096 \times 8\)) times out at \(3{,}600\,\mathrm{s}\) with \(\sim\)​32 MB resident memory and no sign of blowup but no completion either.

The bottleneck on Olivetti is the number of inequalities, not the rank: cdd handles \(400\) inequalities in \(8\)-dim space within \(\sim\)​200 s, while \(4096\) inequalities in the same \(8\)-dim space exhausts \(3{,}600\,\mathrm{s}\). Any successor pipeline that intends to handle the \(m, n \gg r\) regime must either (i) sub-sample the rows passed to DDM, (ii) replace exact extreme-ray enumeration by random sampling, or (iii) abandon the dual representation entirely.

7 Augmented Support and Alt-LP Fallback↩︎

The saturation evidence in Sections 5.2 and 6.1 shows that the obtuseness walk on closed-form \(M_T\) has a structural ceiling. The failing trials have no \(r\)-subset on either side whose \(M_T\) is componentwise nonneg. They are not unreachable by virtue of search budget; they are unreachable by virtue of the uniform-support restriction. We now lift that restriction in the smallest way that still admits a tractable search.

7.1 Augmenting the obtuse cone with \(k = 2\) farthest rays↩︎

For each candidate pair \((T, K)\) visited by the walk, we augment both supports by \(k = 2\) additional ray indices chosen for maximum angular separation from the existing rays of \(T\) (resp. \(K\)). Concretely, for an index \(i \notin T\) we score \[\mathrm{maxcos}(i; T) \;:=\; \max_{j \in T} \frac{|\langle R[:, i], R[:, j] \rangle|}{\|R[:, i]\| \, \|R[:, j]\|}, \label{eq:maxcos}\tag{5}\] and pick the two \(i \notin T\) with the smallest \(\mathrm{maxcos}\) (most isolated from \(T\)). The augmented index set \(T' = T \cup \{i_1, i_2\}\) has size \(r + 2\); we extract the corresponding \(r \times (r + 2)\) basis \(R_{T'} = R_{\mathrm{orth}}[:, T']\). Symmetrically for \(K'\) on the other side.

7.2 Alt-LP for \(\mu, \nu \in \mathbb{R}^{(r+k) \times r}_{\ge 0}\)↩︎

On the augmented bases the witness coefficients become \(\mu, \nu \in \mathbb{R}^{(r+k) \times r}_{\ge 0}\), and the joint factorisation constraint takes the bilinear form \[R_{T'} \, \mu \, \nu^{\!\top} R_{K'}^{\!\top} \;=\; I_r, \label{eq:joint-aug}\tag{6}\] with \(r^2\) scalar equations in \(2(r + k) r\) unknowns. We solve 6 by alternation: fix \(\nu\), solve the linear system for \(\mu\) via slack-LP (Section 2); fix \(\mu\), solve for \(\nu\) similarly. Each subsolve is an LP of size \(r^2\) equations in \((r + k) r\) variables plus \(2 r^2\) slack variables — e.g.\(16\) equations in \(24 + 32 = 56\) variables at \(r = 4, k = 2\). We initialise \(\mu^{(0)} = [\max(M_T, 0); \mathbf{0}]\), \(\nu^{(0)} = [I_r; \mathbf{0}]\) so iteration \(0\) is the original uniform-support solution; the alt-LP only deviates from it as needed.

7.3 Hybrid algorithm↩︎

The full hybrid embeds the alt-LP as a fallback inside the obtuseness walk, used only when the closed-form \(M_T\) is infeasible at the current pair. A per-trial cap \(\mathrm{maxAltLpAttempts} = 1{,}000\) bounds worst-case wall time on instances where many pairs need to be tried.

7.4 Empirical result↩︎

We rerun the \(m = n = 10\), \(r \in \{4, 5, 6\}\), \(100\)-trial Monte Carlo with Algorithm 2 and compare against every prior variant in Table 4.

7.4.0.1 Breakdown.

Within the hybrid, the path that closed each successful trial splits as \(32\) Mt-path \(+ 67\) alt-LP-path at \(r = 4\); \(23 + 72\) at \(r = 5\); \(14 + 61\) at \(r = 6\). The Mt-path count is lower than the Mt-only top-\(200\) baseline (\(79, 85, 58\)) because alt-LP often succeeds at an earlier-ranked pair before the walk reaches the Mt-feasible one later in the order; the alt-LP “steals” wins from Mt. The total (\(99, 95, 75\)) is what matters.

7.4.0.2 Cost.

The hybrid is dominated by the alt-LP fallback at \(r = 6\), where the median number of alt-LP attempts per trial is \(238.5\) (max \(1{,}000\), the cap). Per-trial wall is \(\sim 1\,\mathrm{s}\) at \(r = 4\), \(\sim 41\,\mathrm{s}\) at \(r = 5\), and \(\sim 51\,\mathrm{s}\) at \(r = 6\). Total sweep time was \(\sim 2\,\mathrm{h}\,39\,\mathrm{m}\) on a \(2024\) Apple-silicon laptop. Cone reconstruction error stays at or near machine precision when the hybrid succeeds: medians \(1.3 \times 10^{-15}\), \(5.5 \times 10^{-15}\), \(1.9 \times 10^{-14}\) at \(r = 4, 5, 6\); minima \(\sim 3 \times 10^{-16}\) at every \(r\); maxima \(4 \times 10^{-13}\), \(1.6 \times 10^{-12}\), \(1.2 \times 10^{-11}\) (the alt-LP-path successes can settle right at the \(10^{-10}\) feasibility tolerance rather than to full machine precision).

7.4.0.3 Reading.

The augmented + alt-LP fallback is the first variant in our sequence of experiments that pushes any \(r\) above \(\sim 90\%\) success at this size. The geometric story is simple: a trial in which no \(r\)-subset is near-orthogonal can still have a \((r + 2)\)-subset that admits a feasible \(\mu, \nu\) once we let both sides slacken to \((r + 2)\) atoms per column-bundle. This makes the heuristic substantially more robust, at the cost of introducing an iterative inner LP solve.

7.5 The hybrid at \(m = n = 15\)↩︎

Does the same hybrid recover the failures of the Mt-only \(m = n = 15\) sweep (Table 3)? We ran Algorithm 2 unchanged at \(m = n = 15\) for \(r \in \{4, 5, 6\}\), with the same \(100\)-trial Monte Carlo and \(\mathrm{maxTries} = 200\), \(\mathrm{maxAltLpAttempts} = 1{,}000\) caps used at \(m = n = 10\). Results in Table 5.

The reading is that the hybrid breaks the structural ceiling at the saturation curve’s edge (\(m = n = 10, r = 4, 5, 6\) and \(m = n = 15, r = 4, 5\)), but does not extend the curve into regimes where the underlying cones are uniformly skinny (\(m = n = 15, r \ge 6\)). For the latter, a deeper augmentation (\(k > 2\)) or a fundamentally different relaxation of the uniform-support assumption is needed.

8 Computational Cost Summary↩︎

9 Related Work↩︎

9.0.0.1 Separable / geometric NMF.

The separability assumption, that each cone ray is itself a column of \(A\), makes exact NMF polynomially tractable [5]–[8]. The Successive Projections Algorithm [13] and its descendants exploit this structure. Our toolkit makes no separability assumption; on separable inputs it is overkill, but on the non-separable random rank-\(r\) inputs of Section 5 those algorithms do not return exact witnesses.

9.0.0.2 Algebraic NMF.

The nonneg-rank decision problem is in the first-order theory of real-closed fields [11], [12], [14], [15] and admits Gröbner-basis attacks [16]. The squared-variable substitution we use (\(u_{i,j}^2 \to \alpha_{i,j}\)) is a standard real-algebraic device. The geometric prefilter via Proposition [prop:Mt] is, to our knowledge, original to this work as a way of cutting the combinatorial enumeration cost. The fixed-points of the Mathematica Reduce solver on these systems have been studied geometrically [17], [18].

9.0.0.3 Alternating algorithms.

Multiplicative updates [1], HALS [2], and block coordinate descent [3], [19] are the workhorses of approximate NMF. Our alt-LP differs in two respects: it is applied to an \(r \times (r + k)\) augmented support rather than the full \(m \times r\) matrix, and it minimises slack rather than reconstruction error. The slack-LP plateau noted by [2] for \(\ell_1\)-loss alternating algorithms reappears in our fully-unrestricted alt-LP (Section 10) but is largely avoided in the augmented variant where most of the work is still done by the closed-form \(M_T\).

10 Limitations and Open Problems↩︎

1. DDM scalability. Our pipeline calls cddlib on \(m \times r\) and \(n \times r\) inequality systems. At \(m \sim 400\) and \(r = 8\) this is tractable (\(\sim 200\,\mathrm{s}\)); at \(n = 4096\) and \(r = 8\) it is not (Section 6.2). Replacing the exact DDM by random ray sampling or an incremental partial enumeration is the most pressing item.

2. Augmentation depth \(k\). We set \(k = 2\) throughout. Sweeping \(k \in \{1, 3, 4, 5\}\) at \(r = 6\) to characterise the marginal value of each extra augmentation ray is a clean next experiment. The trade-off is per-pair alt-LP cost (polynomial in \(k\)) against the number of pairs successfully rescued.

3. Initialisation of alt-LP. Our initial \((\mu^{(0)}, \nu^{(0)})\) is the clipped-\(M_T\) embedding of the uniform-support solution. Replacing this by NNLS or a damped Lagrangian inner solve might break the slack-LP plateau on trials where the current alt-LP fails (the \(75\%\) ceiling at \(r = 6\) may not be the geometric ceiling).

4. Lifting the \(r\)-subset structure entirely. A blind per-column LP support discovery (svd_nonuniform.wl) was prototyped: derive per-column supports \(T_j = \mathrm{supp}(\mathrm{lpNN}(R_{\mathrm{orth}}, R_{\mathrm{orth}}[:, j]))\) via \(r\) blind LPs, then iterate full \(k_1 \times r\) and \(k_2 \times r\) \(\mu, \nu\) slack-LPs without support restriction. On the \(r = 4\) sweep this prototype plateaus at residual \(0.35\)–\(1.12\) and rel.err.\(0.054\)–\(0.55\), worse than every baseline. Replacing the slack-LP by NNLS or a damped Lagrangian is the most direct fix.

5. The joint Gröbner basis at \(r = 4\). The companion script svdGb.wl encodes the bilinear system as a Gröbner-basis problem and admits a one-day budget per surviving \((T, K)\) pair after Mt prefiltering. On the \(r = 4\) instance we left running for \(\sim 20\) hours of CPU time, the first surviving pair’s \(\mathrm{GroebnerBasis}+\mathrm{Reduce}\) pipeline has not completed and the memory footprint oscillates between \(0.25\) GB and \(1.0\) GB. A clean characterisation of which \((T, K)\) pairs admit a polynomial-time GB is open.

11 Conclusion↩︎

This paper develops and characterises a cone-ray pipeline for exact NMF on small-\(m, n\) inputs. The closed-form witness \(M_T = R_T^{-1} (R_K^{\!\top})^{-1}\) under uniform support is the cheap test that anchors everything; the obtuseness-ranked walk gives a principled order in which to try \((T, K)\) candidates; and the \(100\)-trial Monte Carlo characterises this walk’s reach precisely. Top-\(200\) is the structural ceiling for closed-form \(M_T\); doubling to top-\(400\) adds \(\le 2\) successes per \(r\), and enlarging \(m\) to \(15\) makes things worse, not better.

The main contribution beyond characterisation is the augmented + alt-LP fallback: at each pair where \(M_T\) is infeasible, augment both supports by \(k = 2\) farthest rays and alt-LP for \(\mu, \nu \in \mathbb{R}^{(r+k) \times r}_{\ge 0}\). This lifts \(m = n = 10\) success counts from \(79/85/58\) to \(\mathbf{99/95/75}\) at \(r = 4, 5, 6\) — the first variant we have found that decisively breaks the structural ceiling. The remaining failures (one at \(r = 4\), five at \(r = 5\), twenty-five at \(r = 6\)) are good targets for the next augmentation depth or for the NNLS/damped-Lagrangian inner solver improvements suggested in Section 10.

Acknowledgments↩︎

The author thanks Anthropic’s Claude (Opus 4.7, 1M context) for help with software drafts, prose, and the Monte-Carlo characterisation pipeline; all errors and design choices are the author’s responsibility. The full executable toolkit, including all reproducible examples and per-trial JSON results, is released alongside this paper as svd.wl, svd_obtuse.wl, svd_obtuse_m10.wl, svd_obtuse_m10_top5.wl, svd_obtuse_m10_top400.wl, svd_obtuse_m10_aug2.wl, svd_obtuse_m10_aug2_top200.wl, svd_obtuse_m15.wl, svd_obtuse_m15_top400.wl, svd_obtuse_m15_aug2_top200.wl, svd_obtuse_olivetti.wl, svdGb.wl, and svd_nonuniform.wl. The complete repository is publicly available at https://github.com/mithilr-sys/cone-nmf-exact.

References↩︎