The Nondecreasing Rank

TL;DR

This paper introduces the concept of nondecreasing (ND) rank for matrices and tensors, explores its properties, and develops algorithms for low ND rank approximations, with applications to real-world datasets.

Contribution

It defines ND rank for tensors, analyzes its theoretical properties, and proposes a hierarchical ALS algorithm for low ND rank tensor approximation.

Findings

ND rank relates to nonnegative rank under certain conditions.

Not all monotonic tensors have finite ND rank.

Low ND rank factorizations reveal meaningful data insights.

Abstract

In this article the notion of the nondecreasing (ND) rank of a matrix or tensor is introduced. A tensor has an ND rank of r if it can be represented as a sum of r outer products of vectors, with each vector satisfying a monotonicity constraint. It is shown that for certain poset orderings finding an ND factorization of rank $r$ is equivalent to finding a nonnegative rank-r factorization of a transformed tensor. However, not every tensor that is monotonic has a finite ND rank. Theory is developed describing the properties of the ND rank, including typical, maximum, and border ND ranks. Highlighted also are the special settings where a matrix or tensor has an ND rank of one or two. As a means of finding low ND rank approximations to a data tensor we introduce a variant of the hierarchical alternating least squares algorithm. Low ND rank factorizations are found and interpreted for two datasets concerning the weight of pigs and a mental health survey during the COVID-19 pandemic.

Figures

Tables

Table 1. Table 1 : Percentages of flies that died after being exposed to three different types of selenium at various concentrations.

	Concentration
	0	100	200	400
Type I	2.0	27.4	26.7	68.0
Type II	1.4	19.6	41.5	40.3
Type III	2.9	24.4	75.0	96.5

Table 2. Table 2 : Percentage of respondents who indicated that they perceived their mental health to be either fair or poor.

Age Groups	Female
	2016	2017	2018	2019	2020	2021	2022
12-17	6.0	7.8	8.5	8.4	12.9	16.5	21.0
18-34	9.0	10.1	12.1	13.1	15.3	17.8	24.2
35-49	7.6	7.9	8.4	8.2	10.8	14.6	16.4
50-64	8.3	8.1	7.9	7.5	9.3	11.6	14.0
65+	5.7	4.7	5.4	5.5	6.0	6.8	9.2
	Male
12-17	2.9	3.9	5.0	3.8	3.8	7.5	8.7
18-34	6.7	6.6	7.6	10.6	11.2	13.6	16.1
35-49	5.2	6.2	6.8	7.3	9.9	11.6	13.5
50-64	7.3	6.4	7.6	6.8	8.7	9.0	11.7
65+	6.3	6.0	4.5	4.9	4.9	6.7	7.9

Table 3. (a) Age group factors.

Age group	${𝐚}^{\left(1\right)}$	${𝐚}^{\left(2\right)}$
12-17	8.93	8.28
18-34	16.77	8.28
35-49	14.23	5.06
50-64	13.18	4.17
65+	8.93	2.53

Table 4. (a) Age group factors.

Age group	${𝐚}^{\left(1\right)}$	${𝐚}^{\left(2\right)}$
12-17	8.93	8.28
18-34	16.77	8.28
35-49	14.23	5.06
50-64	13.18	4.17
65+	8.93	2.53

Table 5. (b) Year factors.

Year	${𝐛}^{\left(1\right)}$	${𝐛}^{\left(2\right)}$
2016	0.76	0.57
2017	0.75	0.70
2018	0.82	0.78
2019	0.90	0.75
2020	1.02	1.07
2021	1.25	1.34
2022	1.49	1.78

Table 6. (c) Sex factors.

Sex	${𝐜}^{\left(1\right)}$	${𝐜}^{\left(2\right)}$
Female	0.38	1
Male	0.62	0

Equations

$$T_{i_1\dots i_k}=v_{i_1}^{(1)}\cdots v_{i_k}^{(k)},\forall (i_1,\dots ,i_k)\in {\times }_{j=1}^k[p_j],$$

$$\text{rank}(\mathbf{T})=\min \{r:\mathbf{T}=\sum _{i=1}^r{\otimes }_{j=1}^k{\mathbf{v}}^{(ij)}\},$$

$${\text{rank}}_{+}(\mathbf{T})=\min \{r:\mathbf{T}=\sum _{i=1}^r{\otimes }_{j=1}^k{\mathbf{v}}^{(ij)},v_l^{(ij)}\ge 0\forall i,j,l\}.$$

$$\mathcal{C}(\mathcal{P})={\cap }_{x⪯y}\{\mathbf{f}:f_x\le f_y\}\cap \{\mathbf{f}:f_x\ge 0,x\in \mathcal{P}\}.$$

$$\mathrm{NDrank}(\mathbf{T})=\min \{r:\mathbf{T}=\sum _{i=1}^r{\otimes }_{j=1}^k{\mathbf{v}}^{(ij)},{\mathbf{v}}^{(ij)}\in \mathcal{C}({\mathcal{P}}_j),\forall i,j\}.$$

$${\otimes }_{j=1}^k{\mathcal{C}}_j:-\{\mathbf{T}:\exists r,\mathbf{T}=\sum _{i=1}^r{\otimes }_{j=1}^k{\mathbf{v}}^{(ij)},{\mathbf{v}}^{(ij)}\in {\mathcal{C}}_j\}.$$

$$\left[\begin{array}{ccc}1& 1& 1\\ 1& 1& 1\end{array}\right],\left[\begin{array}{ccc}0& 1& 1\\ 0& 1& 1\end{array}\right],\left[\begin{array}{ccc}0& 0& 0\\ 1& 1& 1\end{array}\right],\left[\begin{array}{ccc}0& 0& 0\\ 0& 1& 1\end{array}\right],\left[\begin{array}{ccc}0& 0& 1\\ 0& 0& 1\end{array}\right],\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\end{array}\right],$$

$$\left[\begin{array}{ccc}0& 1& 1\\ 1& 1& 1\end{array}\right],\left[\begin{array}{ccc}0& 0& 1\\ 1& 1& 1\end{array}\right],\left[\begin{array}{ccc}0& 0& 1\\ 0& 1& 1\end{array}\right].$$

$$\{(y_1,y_2,x_3,\dots ,x_k):a_1⪯y_1,b_2⪯y_2\}\cup \{(y_1,y_2,x_3,\dots ,x_k):b_1⪯y_1,a_2⪯y_2\}$$

$$-t_{11}\le 0,-t_{12}\le 0,-t_{21}\le 0,-t_{22}\le 0,$$

$$t_{21}-t_{23}\le 0,t_{12}-t_{32}\le 0,t_{22}-t_{23}\le 0,t_{11}-t_{31}\le 0,$$

$$t_{12}-t_{13}\le 0,t_{22}-t_{32}\le 0,t_{21}-t_{31}\le 0,t_{11}-t_{13}\le 0,$$

$$-t_{11}+t_{13}+t_{31}-t_{33}\le 0,-t_{22}+t_{23}+t_{32}-t_{33}\le 0,$$

$$-t_{21}+t_{23}+t_{31}-t_{33}\le 0,-x_2+t_{13}+t_{32}-t_{33}\le 0,$$

$$t_{11}+t_{12}-t_{13}+t_{21}-t_{22}-t_{31}\le 0,t_{11}+t_{12}-t_{13}-t_{21}+t_{22}-t_{32}\le 0,$$

$$t_{11}-t_{12}+t_{21}+t_{22}-t_{23}-t_{31}\le 0,-t_{11}+t_{12}+t_{21}+t_{22}-t_{23}-t_{32}\le 0,$$

$$-t_{11}-t_{12}+t_{13}-t_{21}+t_{22}+t_{31}-t_{33}\le 0,-t_{11}-t_{12}+t_{13}+t_{21}-t_{22}+t_{32}-t_{33}\le 0,$$

$$-t_{11}+t_{12}-t_{21}-t_{22}+t_{23}+t_{31}-t_{33}\le 0,t_{11}-t_{12}-t_{21}-t_{22}+t_{23}+t_{32}-t_{33}\le 0.$$

$$\left[\begin{array}{ccc}1& -1& 0\\ 1& 1& -1\\ -1& 0& 1\end{array}\right],\left[\begin{array}{ccc}-1& 1& 0\\ 1& 1& -1\\ 0& -1& 1\end{array}\right].$$

$$[({\mathbf{A}}^{(1)}\otimes \cdots \otimes {\mathbf{A}}^{(k)})(\mathbf{T}){]}_{j_1\dots j_k}=\sum _{i_1=1}^{p_1}\cdots \sum _{i_k=1}^{p_k}{A}_{j_1{i}_1}^{(1)}\cdots A_{j_k{i}_k}^{(k)}{T}_{i_1\dots i_k}.$$

$${\mathbf{M}}_{\mathcal{P}}({\mathbf{e}}_x)=\sum _{x⪯y}{\mathbf{e}}_y,\forall \mathbf{x}\in \mathcal{P},$$

$${\mathbf{M}}_{\mathcal{P}}{\mathbf{e}}_{(x_1,\dots ,x_k)}=\sum _{(x_1,\dots ,x_k)⪯(y_1,\dots ,y_k)}{\mathbf{e}}_{(x_1,\dots ,x_k)}=\sum _{x_1⪯y_1}\cdots \sum _{x_k⪯y_k}{\mathbf{e}}_{(y_1,\dots ,y_k)}=({\otimes }_{j=1}^k{\mathbf{M}}_{{\mathcal{P}}_j})({\mathbf{e}}_{(x_1,\dots ,x_k)}),$$

$$[{\mathbf{M}}_{\mathcal{P}}^{-1}{]}_{xy}=\{\begin{array}{ll}1,& x=y\\ -1,& y⋖x\\ 0,& \mathrm{otherwise}\end{array}.$$

$${\mathbf{M}}_{\mathcal{P}}^{-1}=\left[\begin{array}{cccccc}1& 0& 0& 0& \cdots & 0\\ -1& 1& 0& 0& \cdots & 0\\ 0& -1& 1& 0& \cdots & 0\\ ⋮& \ddots & \ddots & \ddots & \ddots & ⋮\\ ⋮& \ddots & \ddots & \ddots & \ddots & 0\\ 0& \cdots & \cdots & 0& -1& 1\end{array}\right],$$

$$[{\mathbf{M}}_{\mathcal{P}}{\mathbf{M}}_{{\mathcal{P}}^{-1}}{]}_{\cdot y}$$

$$=\sum _{x\in [y,\infty )}{\mathbf{e}}_x-\sum _{z:y⋖z}\sum _{x\in [z,\infty )}{\mathbf{e}}_x={\mathbf{e}}_y.$$

$$[({\otimes }_{j=1}^k{\mathbf{M}}_{{\mathcal{P}}_j}^{-1})(\mathbf{T}){]}_{i_1\dots i_k}=\sum _{j_1=0}^1\cdots \sum _{j_k=0}^1(-1{)}^{{\sum }_{l=1}^k{j}_l}{T}_{i_1-j_1,\dots ,i_k-j_k},$$

$$({\otimes }_{j=1}^2{\mathbf{M}}_{{\mathcal{P}}_j}^{-1})(\mathbf{T})$$

$$=\left[\begin{array}{ccc}{t}_{11}& t_{12}-t_{11}& t_{13}-t_{12}\\ t_{21}-t_{11}& t_{22}-t_{21}-t_{12}+t_{11}& t_{23}-t_{22}-t_{13}+t_{12}\end{array}\right].$$

$$\sum _{i_1\in \{x_1,y_1\}}\cdots \sum _{i_k\in \{x_k,y_k\}}(-1{)}^{{\sum }_{j=1}^k\delta (i_j)}{T}_{i_1\dots i_k}\ge 0,$$

$$\{\mathbf{T}:{\otimes }_{j=1}^k({\mathbf{e}}_{y_j}^{\ast }-{\mathbf{e}}_{x_j}^{\ast })(\mathbf{T})=\sum _{i_1\in \{x_1,y_1\}}\cdots \sum _{i_k\in \{x_k,y_k\}}(-1{)}^{{\sum }_{j=1}^k\delta (i_j)}{T}_{i_1\dots i_k}\ge 0\}.$$

$$\text{maxNDrank}=\underset{r<\infty }{\sup }\{{\mathcal{N}}_r\ne \varnothing \}.$$

$$\sum _{i=1}^r{\otimes }_{j=1}^k{\mathbf{v}}^{(ij)}$$

$$=\sum _{l_1=1}^{q_1}\cdots \sum _{l_{k-1}=1}^{q_{k-1}}(\underset{j=1}{\overset{k-1}{⨂}}{\mathbf{w}}^{(l_{j}j)}\otimes (\sum _{i=1}^r\sum _{l_k=1}^{q_k}\left(\prod _{j=1}^k{\alpha }_{ij{l}_j}\right){\mathbf{w}}^{(l_{k}k)}\left)\right),$$

$$F\backslash (\bigcup _{J:|J|<p-1,\forall j\in J{\mathbf{v}}_j\in F}\text{cone}({\mathbf{v}}_j:j\in J)),$$

$$0={\mathbf{h}}^{⊺}\sum _{i\in I}{\alpha }_i({\mathbf{v}}_i\otimes {\mathbf{v}}_i){\mathbf{h}}^{\prime }={\mathbf{h}}^{⊺}(\sum _{i=1}^r\sum _{j=1}^{2^{p-1}}\sum _{l=1}^{2^{p-1}}{\lambda }_{ij}{\mu }_{il}({\mathbf{v}}_j\otimes {\mathbf{v}}_l)){\mathbf{h}}^{\prime }\ge \sum _{i=1}^r{\lambda }_{ij}{\mu }_{il}({\mathbf{h}}^{⊺}{\mathbf{v}}_j)({\mathbf{v}}_l^{⊺}{\mathbf{h}}^{\prime }).$$

Full text

Keywords: contingency table; hierarchical alternating least squares; Möbius inversion; nonnegative matrix factorization; partially ordered set; polytope; simplicial cone; tensor; rank.

1 Introduction

Matrix and tensor factorizations are indispensable tools in data science, both for finding easily interpretable representations of complex data, and for obtaining more stable estimates of estimands of interest. An important example of a matrix factorization technique is nonnegative matrix factorization (NMF) [29], which has been leveraged in many applications ranging from hyperspectral imaging [58] to recommender systems [39]. In a nonnegative matrix factorization a matrix is decomposed into a product of two, low-rank, nonnegative matrices. The rows and columns of the matrix factors often have intuitive interpretations as atomic units that the dataset is built from.

Nonnegative matrix factorizations have a statistical motivation as they also arise in the study of contingency tables as mixture models of independent random variables [22, Ch 4]. One strand of literature within the field of algebraic statistics has been to understand the semialgebraic geometry of such mixture models. Despite being simple to formulate, the geometry of NMF is involved, even in the rank-$3$ case [23, 33].

The present paper is concerned with a generalization of NMF; in addition to nonnegativity constraints, we seek matrix and tensor factorizations that fulfill monotonicity constraints with respect to a user-specified partial order. Termed a nondecreasing (ND) factorization, much of the theory for nonnegative factorizations will hold in this generalized setting. There are some notable exceptions where the theory for ND factorizations differs. For instance, the maximum ND rank of a matrix can be significantly larger than the maximum nonnegative rank.

Monotonicity constraints with respect to a partially ordered set feature prominently in order constrained statistical inference [49]. Order constraints often arise in hypothesis testing problems where the efficacy of a collection of treatments are compared to determine which treatment appears to produce the best results [50]. When observing functional data order constraints are also relevant. It may for example be natural to require that a function be monotone increasing [44]. Monotonicity constraints for matrix or tensor data have been previously modeled by multi-way, order constrained, analyses of variance (ANOVA) [28]. ND factorization, while similar, is not equivalent to order constrained ANOVA.

An overview of this paper is as follows: Section 3 provides background material and Section 4 provides further introduction and motivation for the concept of the ND rank. Sections 5 and 6 address the questions of when exact ND factorizations exist and when they can be reframed as nonnegative factorizations up to a linear change in coordinates. The geometry of order cones is also explicated here. Section 7 provides exact expressions and bounds for the maximum ND rank in certain cases. It is shown that knowing the maximum ND rank is enough to fully determine the typical ND ranks. Section 8 shows that if the rank of a matrix is one or two and it has a finite ND rank, then the ND rank is also equal to one or two respectively. Section 9 shows that the border ND rank is equivalent to the ND rank, a fact which is relevant for optimization. Section 10 proposes an ND hierarchical least squares optimization algorithm. The paper concludes with an application section illustrating how ND factorizations are fruitful for uncovering structure in data.

2 Notation

Some of the notation used in this work is listed here for easy reference. Partially ordered sets (posets) are represented by $\mathcal{P}$, the order cone associated with this poset is denoted by $\mathcal{C}(\mathcal{P})$, and the set of tensors with ND rank at most $r$ is $\mathcal{N}_{\leq r}$. The dual of a cone $\mathcal{C}$ is $\mathcal{C}^{*}$. A curly inequality $\prec$ represents inequality with respect to a poset while $x\lessdot y$, means that the element $y$ covers $x$ in the poset. $\mathbb{R}^{p}_{+}$ is equal to all vectors in $\mathbb{R}^{p}$ with nonnegative entries and ${\mathbf{e}}_{i}$ is the $i$th standard basis vector. Vectors, matrices, and tensors are all displayed in bold font, while scalars are not. Given any finite-dimensional vector spaces $V_{1},\ldots,V_{k}$ and sets $S_{i}\subseteq V_{i}$, $\otimes_{i=1}^{k}V_{i}$ is the tensor product space of the $V_{i}$ and $\otimes_{i=1}^{k}S_{i}$ is the subset $\cup_{r=1}^{\infty}\{\sum_{i=1}^{r}\otimes_{j=1}^{k}{\mathbf{s}}^{(ij)}:{\mathbf{s}}^{(ij)}\in S_{j}\}$ of this tensor product space. We will use $k$ to refer to the number of factors in the tensor product and take $p_{i}\coloneq\text{dim}(V_{i})$, where usually $V_{i}=\mathbb{R}^{\mathcal{P}_{j}}$ with $\mathcal{P}_{j}$ a poset containing $p_{j}$ elements. Shorthand for the set $\{1,\ldots,p\}$ is $[p]$. The standard inner product between two order-$k$ tensors is $\langle{\mathbf{S}},{\mathbf{T}}\rangle=\sum_{i_{1}=1}^{p_{1}}\cdots\sum_{i_{k}=1}^{p_{k}}S_{i_{1}\ldots i_{k}}T_{i_{1}\ldots i_{k}}$. This induces the Frobenius norm $\|{\mathbf{T}}\|_{F}^{2}=\langle{\mathbf{T}},{\mathbf{T}}\rangle$.

3 Preliminaries on Low-rank Tensors and Convex Geometry

For the purposes of this work a tensor ${\mathbf{T}}\in\mathbb{R}^{p_{1}\times\cdots\times p_{k}}$ will be viewed as an array of real numbers $T_{i_{1}\ldots i_{k}}$ where $i_{j}$ ranges from $1$ to $p_{j}$ for all $j$. Equivalently, a tensor ${\mathbf{T}}$ can be viewed as a function on the index set $\times_{j=1}^{k}[p_{j}]$ with ${\mathbf{T}}(i_{1},\ldots,i_{k})\coloneq T_{i_{1}\ldots i_{k}}$. We say that such a tensor has order $k$ and dimension $(p_{1},\ldots,p_{k})$. Order one tensors are vectors and order two tensors are matrices. A large portion of the applications of this work are to matrix-valued data and a reader unacquainted with tensors may restrict their focus to matrices. Modes of the tensor refer to one of the $k$ index sets $[p_{j}]$, and a fibre of this tensor is a vector that holds every index of ${\mathbf{T}}$ fixed except for one. For example, ${\mathbf{S}}_{\bullet 4}\in\mathbb{R}^{p_{1}}$ represents the fourth column of a matrix — a mode-one fibre, while ${\mathbf{T}}_{3\,2\bullet 4}\in\mathbb{R}^{p_{3}}$ is a mode-three fibre of an order four tensor.

A special class of tensors are tensors that have rank-one, meaning that they can be written multiplicatively as

[TABLE]

where every $({\mathbf{v}}^{(j)})^{\intercal}=(v_{1}^{(j)},\ldots,v_{p_{j}}^{(j)})$ is a vector in $\mathbb{R}^{p_{j}}$. The tensor product notation $\otimes$ is used to denote the above rank-one structure: ${\mathbf{T}}=\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}$. The real (CP [32, Sec 3]) rank of a tensor is defined as

[TABLE]

where every ${\mathbf{v}}^{(ij)}\in\mathbb{R}^{p_{j}}$. That is, it is the smallest number of rank-one tensors required to express ${\mathbf{T}}$ as a sum of rank-one tensors. Whenever rank is referred to in this work it is taken to mean rank over the real numbers. The nonnegative rank [29] places extra nonnegativity constraints on the vectors ${\mathbf{v}}^{(ij)}$ and is defined as

[TABLE]

In general, the nonnegative rank of a tensor may exceed its rank [17, 55].

We will be examining low-rank decompositions where the vectors appearing in the rank-one factors satisfy a conical constraint. A brief summary of some standard ideas from convex geometry are provided below, where the reader is referred to [12, 59] for more details.

A (convex) cone $\mathcal{C}\subset\mathbb{R}^{p}$ satisfies the property that for all ${\mathbf{x}},{\mathbf{y}}\in\mathcal{C}$ and $\alpha,\beta\geq 0$ the vector $\alpha{\mathbf{x}}+\beta{\mathbf{y}}$ is in $\mathcal{C}$. The conical hull of the vectors ${\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{m}$ is defined as $\{\sum_{i=1}^{m}\alpha_{i}{\mathbf{v}}_{i}:\alpha_{i}\geq 0\}$. Any cone that can be expressed in this way is called polyhedral with the representation referred to as a vertex or $\mathcal{V}$-representation. Polyhedral cones can also equivalently be expressed in a halfspace or $\mathcal{H}$-representation, where a cone is defined by finitely many halfspace constraints: $\{{\mathbf{x}}:{\mathbf{a}}_{i}^{\intercal}{\mathbf{x}}\geq 0,i=1,\ldots,l\}$. The dual cone $\mathcal{C}^{*}$ is closely related to the $\mathcal{H}$-representation as it is defined as the closed, convex set $\mathcal{C}^{*}=\{{\mathbf{a}}:{\mathbf{a}}^{\intercal}{\mathbf{x}}\geq 0,\forall{\mathbf{x}}\in\mathcal{C}\}$. A cone is said to be pointed if $\mathcal{C}\cap(-\mathcal{C})=\{{\mathbf{0}}\}$, solid if it has a non-empty interior, and proper if it is pointed, solid, and closed. Polyhedral cones are pointed if and only if the ${\mathbf{a}}_{i}$ appearing in the $\mathcal{H}$-representation span $\mathbb{R}^{p}$, are solid if and only if the ${\mathbf{v}}_{i}$ appearing in the $\mathcal{V}$-representation span $\mathbb{R}^{p}$, and are always closed. Whenever $\mathcal{C}$ is polyhedral or proper the dual cone $\mathcal{C}^{*}$ inherits these properties respectively.

One of the simplest types of polyhedral cones are simplicial cones that are generated by $p$ linearly independent vectors ${\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{p}\in\mathbb{R}^{p}$; other classes of proper cones are generated by more than $p$ vectors. The nonnegative orthant $\mathbb{R}^{p}_{+}=\{\sum_{i=1}^{p}\alpha_{i}{\mathbf{e}}_{i}:\alpha_{i}\geq 0\}$ is simplicial as it is generated by the standard basis vectors ${\mathbf{e}}_{i}$. A face $\mathcal{F}\subset\mathcal{C}$ of a cone is a subset where if ${\mathbf{x}}+{\mathbf{y}}\in\mathcal{F}$ and ${\mathbf{x}},{\mathbf{y}}\in\mathcal{C}$ then ${\mathbf{x}},{\mathbf{y}}\in\mathcal{F}$. One-dimensional faces, referred to as extremal rays, have the form $\{\alpha{\mathbf{v}}:\alpha\geq 0\}$ of a ray protruding from ${\mathbf{0}}$ in the direction ${\mathbf{v}}$. For brevity we will at times refer to ${\mathbf{v}}$ as an extremal ray or as being extremal. The extremal rays of $\mathbb{R}_{+}^{p}$ are exactly the coordinate axes $\{\alpha{\mathbf{e}}_{i}:\alpha\geq 0\}$. Any face of $\mathbb{R}^{p}_{+}$ is equal to $\{\sum_{i\in I}\alpha_{i}{\mathbf{e}}_{i}:\alpha_{i}\geq 0\}=\{{\mathbf{x}}\in\mathbb{R}^{p}_{+}:x_{i}=0,\forall i\notin I\}$ for some set of indices $I\subset[p]$. The set of faces of a polyhedral cone forms a partially ordered set under set inclusion. For example, the face $\{\alpha_{1}{\mathbf{e}}_{1}+\alpha_{2}{\mathbf{e}}_{2}:\alpha_{1},\alpha_{2}\geq 0\}$ contains the two extremal rays $\{\alpha{\mathbf{e}}_{1}:\alpha\geq 0\},\{\alpha{\mathbf{e}}_{2}:\alpha\geq 0\}$. As this example illustrates, any face in a proper cone is equal to the conical hull of the extremal rays contained in the face. In particular, if $\mathcal{C}=\text{cone}({\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{m})$ then the faces of $\mathcal{C}$ are equal to $\text{cone}({\mathbf{v}}_{i_{j}}:i_{j}\in I)$ for some index set $I$, but not every index set $I\subset[m]$ gives rise to a face. A facet of a cone is a face that has dimension equal to the dimension of the cone minus one.

4 The Nondecreasing Rank

Going beyond nonnegativity constraints, in many applications a tensor that arises from data may also possess some kind of ordering constraint. Consider the following dose-response data in Table 1 that displays the percentage of flies that died in response to exposure to different types of selenium at different concentrations [31, 46]. As might be expected, there is a general increasing trend within each row of the table, as more flies expire when selenium is present in higher concentrations. It is also apparent from the table that the third type of selenium appears more harmful than the first two types, while the first type may be more harmful than the second, but this is much less clear.

If a concise representation of this data was sought, a simple assumption would be that the dose-response functions $f_{i}:\{0,100,200,400\}\rightarrow[0,100]$ for each type of selenium satisfy $f^{(i)}=c_{i}f$ for some nondecreasing function $f$ and scalars $c_{3}\geq c_{1},c_{2}\geq 0$, with the latter assumption coming from the observation that the third type of selenium appears to be most harmful. This representation amounts to assuming that the dose-response matrix ${\mathbf{T}}\approx{\mathbf{c}}{\mathbf{f}}^{\intercal}$ is rank-one and that the nonnegative vectors ${\mathbf{c}}$ and ${\mathbf{v}}$ are congruous with the respective orderings $\text{Type I, Type II}\leq\text{Type III}$, and $0\leq 100\leq 200\leq 400$ that are imposed on the row and column variables of the matrix. If ${\mathbf{c}}$ is scaled so that $\tfrac{1}{3}(c_{1}+c_{2}+c_{3})=1$ then the vector ${\mathbf{f}}$ could be interpreted as the average dose-response curve across all three types of selenium, while ${\mathbf{c}}$ reflects the multiplicative, noxious effects of the different varieties of selenium.

Formally, in the selenium example we have associated two partially ordered sets (posets) with the rows and columns of the table, namely $\mathcal{P}_{1}=\{\mathrm{Type\,I,Type\,II,Type\,III}\}$ and $\mathcal{P}_{2}=\{0,100,200,400\}$. The vectors ${\mathbf{c}}$ and ${\mathbf{f}}$ can be viewed as functions defined on $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ respectively. To say that both ${\mathbf{c}}$ and ${\mathbf{f}}$ respect the given partial orderings is equivalent to saying that these vectors lie in the order cones $\mathcal{C}(\mathcal{P}_{1})$ and $\mathcal{C}(\mathcal{P}_{2})$:

Definition 1** (Order cone and order polytope of functions).**

Let $\mathcal{P}$ be a finite poset and take ${\mathbf{f}}\in\mathbb{R}^{\mathcal{P}}$ be a real-valued function defined on this poset. The function ${\mathbf{f}}$ is said to be nondecreasing if $f_{x}\leq f_{y}$ whenever $x\preceq y$ in $\mathcal{P}$. The order cone $\mathcal{C}(\mathcal{P})$ consists of all nonnegative, nondecreasing functions on $\mathcal{P}$, while the order polytope $\mathcal{O}(\mathcal{P})$ consists of all tensors in $\mathcal{C}(\mathcal{P})$ with $f_{x}\leq 1$ for all $x\in\mathcal{P}$ [51].

An order cone is a proper, polyhedral cone that has the $\mathcal{H}$-representation

[TABLE]

More will be said about the cone structure of $\mathcal{C}(\mathcal{P})$ in the next section.

It may be the case that a rank-one factorization of the form ${\mathbf{c}}{\mathbf{f}}^{\intercal}$ is not sufficient to capture all of the structure inherent to the selenium dose-response matrix; more general low-rank representations are sought. Pairing a monotonicity requirement with a low-rank factorization we now introduce the central notion of interest in this work:

Definition 2** (Nondecreasing (ND) rank).**

Let ${\mathbf{T}}\in\mathbb{R}^{\mathcal{P}_{1}\times\cdots\times\mathcal{P}_{k}}$ be a tensor defined on the Cartesian product index set $\mathcal{P}\coloneq\mathcal{P}_{1}\times\cdots\times\mathcal{P}_{k}$, where each $\mathcal{P}_{j}$ is a finite poset. Any decomposition of the form ${\mathbf{T}}=\sum_{i=1}^{r}\otimes_{j=1}^{k}{\mathbf{v}}^{(ij)}$ with ${\mathbf{v}}^{(ij)}\in\mathcal{C}(\mathcal{P}_{j})$ is referred to as a nondecreasing decomposition or factorization. The nondecreasing rank of a tensor ${\mathbf{T}}$ is the minimum number of rank-one factors needed in a nondecreasing factorization:

[TABLE]

If ${\mathbf{T}}$ does not possess a nondecreasing factorization then it has an infinite ND rank. The set of all tensors with finite ND rank, ND rank equal to $r$, and ND rank less than or equal to $r$ are respectively denoted by $\mathcal{N}_{<\infty}(\mathcal{P})$, $\mathcal{N}_{r}(\mathcal{P})$, and $\mathcal{N}_{\leq r}(\mathcal{P})$. The poset $\mathcal{P}$ will generally be dropped from this notation.

The ND rank refines the nonnegative rank (2) by possibly adding extra constraints, implying that $\mathrm{rank}_{+}({\mathbf{T}})\leq\mathrm{NDrank({\mathbf{T}}})$. We note that the terminology of monotone rank has appeared in the literature [2], although it is synonymous with the nonnegative rank over an ordered field. Other very general notions of rank that encompass the ND rank as a special case are that of the X-rank [43] and the atomic cone rank [24].

The tensor ${\mathbf{T}}$ appearing in the definition of ND rank can be viewed as a function on the Cartesian product of posets $\mathcal{P}\coloneq\times_{j=1}^{k}\mathcal{P}_{j}$. The set $\mathcal{P}$ can itself be given the structure of a poset, called the product poset, by defining the ordering with $(x_{1},\ldots,x_{k})\preceq(y_{1},\ldots,y_{k})$ if and only if $x_{j}\preceq y_{j}$ in $\mathcal{P}_{j}$ for every $j$. In the selenium example we have that $(\mathrm{Type\,II},300)\preceq(\mathrm{Type\,III},400)$ but $(\mathrm{Type\,III},0)\npreceq(\mathrm{Type\,I},400)$. A helpful way to visualize posets is by their Hasse diagrams. A Hasse diagram is a directed graph with nodes given by the elements of the poset and arrows drawn from an element $x$ to an element $y$ if and only if $y$ covers $x$, meaning that $x\prec y$ and there does not exist an element $z$ with $x\prec z\prec y$. The covering relation is denoted by $x\lessdot y$. Figures 1 and 2 illustrate the Hasse diagrams of the partial orderings $\mathcal{P}_{1}$, $\mathcal{P}_{2}$, and $\mathcal{P}_{1}\times\mathcal{P}_{2}$ in the selenium example. We see that there is an arrow from [math] to $100$ as $0\lessdot 100$, but not from [math] to $200$ since $100$ lies in between [math] and $200$. The Hasse diagram of the poset $\mathcal{P}_{2}$ is a path and has the additional property of being totally ordered. Totally ordered posets are referred to as chains: for any two elements $x,y\in\mathcal{P}_{2}$ either $x\preceq y$ or $y\preceq x$. The poset $\mathcal{P}_{1}$ is not a chain since Type I is neither less than or greater than Type II.

It is worth contemplating why there is a requirement that vectors must be nonnegative in the definition of the order cone, and accordingly in ND factorizations. One motivation for this is that if ${\mathbf{T}}=\sum_{i=1}^{r}\otimes_{j=1}^{k}{\mathbf{v}}^{(ij)}$ is an ND factorization then ${\mathbf{T}}\in\mathcal{C}(\mathcal{P})$ is guaranteed to also respect the order of the product poset $\mathcal{P}$. If we contrarily permitted a nondecreasing ${\mathbf{v}}^{(ij)}$ to contain negative entries the tensor $\otimes_{j=1}^{k}{\mathbf{v}}^{(ij)}$ would not necessarily be nondecreasing over $\mathcal{P}$.

5 Geometry of the Order Cone and the Existence of Nondecreasing Factorizations

Any nonnegative tensor has a nonnegative factorization. However, the situation differs for monotone tensors; if ${\mathbf{T}}\in\mathcal{C}(\mathcal{P})$ it is not necessarily the case that ${\mathbf{T}}$ has a nondecreasing factorization. In this section and the next we highlight conditions that ensure that a tensor has a finite ND rank. The geometry of order cones and the finite-rank cone will also be examined.

Definition 3** (Projective Tensor Product of Cones).**

Given convex cones $\mathcal{C}_{j}$, $j=1,\ldots,k$, the projective tensor product [41] of the $\mathcal{C}_{j}$s is equal to

[TABLE]

When $\mathcal{C}_{j}=\mathcal{C}(\mathcal{P}_{j})$ in the above definition the tensor product of cones is the same as the set $\mathcal{N}_{<\infty}(\times_{i=1}^{k}\mathcal{P}_{j})$. The projective tensor product is a proper, polyhedral cone of dimension $\prod_{j=1}^{k}p_{j}$ whenever every $\mathcal{C}_{j}\subset\mathbb{R}^{p_{j}}$ is proper and polyhedral [41]. As the projective tensor product is equal to $\text{cone}(\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}:{\mathbf{v}}^{(j)}\in\mathcal{C}_{j})$ the extremal rays must be rank-one tensors. The following result [18, Thm 3.22] identifies all such extremal rays.

Lemma 1**.**

The extremal rays of $\otimes_{j=1}^{k}\mathcal{C}_{j}$ consist exactly of the rank-one tensors $\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}$ where every ${\mathbf{v}}^{(j)}\in\mathcal{C}_{j}$ is extremal.

In general, the facial structure of $\otimes_{j=1}^{k}\mathcal{C}_{j}$ is more complicated. An exception is in the special case where almost every $\mathcal{C}_{j}$ is simplicial and the dual cone can be described [6, 5, 18]:

Theorem 1**.**

If at most one of the $\mathcal{C}_{j}$ cones is not simplicial then $\big(\otimes_{j=1}^{k}\mathcal{C}_{j}\big)^{*}=\otimes_{j=1}^{k}\mathcal{C}_{j}^{*}$. If every $\mathcal{C}_{j}$ is also polyhedral then the facets of $\otimes_{j=1}^{k}\mathcal{C}_{j}$ have the form $\{{\mathbf{T}}:\langle{\mathbf{T}},\otimes_{j=1}^{k}{\mathbf{h}}^{(j)}\rangle=0\}\cap\otimes_{j=1}^{k}\mathcal{C}_{j}$ where every ${\mathbf{h}}^{(j)}$ is extremal in $\mathcal{C}_{j}^{*}$.

Using Lemma 1, a $\mathcal{V}$-representation of $\mathcal{N}_{<\infty}$ can be obtained from the extremal rays of each $\mathcal{C}(\mathcal{P}_{j})$. The next result [26, Thm 7] describes the extremal rays of any order cone in terms of upsets. An upset $\mathcal{U}\subseteq\mathcal{P}$ is a subset with the property that if $x\in\mathcal{U}$ and $x\preceq y$ then $y\in\mathcal{U}$. An upset is connected if the subgraph of the Hasse diagram corresponding to the upset is connected.

Lemma 2** ([26]).**

The extremal rays of $\mathcal{C}(\mathcal{P})$ consist vectors of the form $\sum_{i\in I}{\mathbf{e}}_{i}$, where $I$ is a non-empty, connected upset of $\mathcal{P}$.

Order cones that are simplicial will be of special interest due to Theorem 1 and the linear equivalence of simplicial cones and the nonnegative orthant. To describe orders that have simplicial order cones we first introduce some terminology borrowed from the literature on directed graphical models [42].

Definition 4**.**

A collider in the Hasse diagram of a poset $\mathcal{P}$ refers to a subgraph of the form

cab

Equivalently, a collider is a set of three elements $a,b,c\in\mathcal{P}$ where $a\lessdot c$ and $b\lessdot c$.

Theorem 2**.**

The cone $\mathcal{C}(\mathcal{P})$ is simplicial if and only if the Hasse diagram does not contain any colliders. A Hasse diagram has no colliders if and only if as an undirected graph it is a disjoint union of trees where each tree has edges directed away from one its nodes.

Proof.

If $\mathcal{P}$ has $p$ elements, by Lemma 2, the cone $\mathcal{C}(\mathcal{P})$ is simplicial if and only if the number of connected upsets of $\mathcal{P}$ is $p$. For every $x\in\mathcal{P}$, the intervals $[x,\infty)=\{y:x\preceq y\}$ are distinct connected upsets. Thus, the number of connected upsets is at least $p$. To show that the order cone is simplicial it suffices to show that the number of connected upsets within each connected component of size $m$ in the Hasse diagram is equal to $m$.

Assume that the Hasse diagram contains no colliders. Any upset $\mathcal{U}$ has the form $\mathcal{U}=\cup_{u\in\mathcal{U}}[u,\infty)$. If $\mathcal{U}$ contains two or more minimal elements $u_{1},u_{2}$ the upset $\mathcal{U}$ will not be connected since the only undirected path in $\mathcal{U}$ from $u_{1}$ to $u_{2}$ must be directed away from both $u_{1}$ and $u_{2}$ since these are minimal elements in $\mathcal{U}$. Any such path must include a collider, so no such path exists. Every connected upset $\mathcal{U}$ must therefore have a minimum element $u$, where $\mathcal{U}=[u,\infty)$. The number of connected upsets in a connected component of size $m$ in this Hasse diagram is equal to $m$.

Conversely, assume that $\mathcal{P}$ has a connected component of the Hasse diagram of size $m$ that has a collider. There exists three elements $a,b,c$ with $a,b$ incomparable and $a\prec c$, $b\prec c$. The upset $a\cup b\cup[c,\infty)$ is connected and does not equal $[x,\infty)$ for any $x\in\mathcal{P}$. This shows that this connected component of $\mathcal{P}$ has at least $m+1$ connected upsets.

The last statement that the absence of colliders implies the given tree structure is proven by noting that if the undirected graph associated with the Hasse diagram contains a cycle, the only way to orient the edges in a cycle to avoid the appearance of a collider is to create a directed cycle. This contradicts the fact that Hasse diagrams are directed acyclic graphs. If a tree does not contain a minimum element that all edges are directed away from there exists at least two distinct minimal elements $a,b$ in the tree. As the tree is connected there exists a path between $a$ and $b$ with arrows directed away from both $a$ and $b$. There is no way to orient the edges of such a path without obtaining a collider. The converse direction that unions of trees with the given structure contain no colliders is apparent. ∎

Figure 3 illustrates the tree structure described in Theorem 2, while Figure 4 (a) and (b) respectively illustrate order cones that are and are not simplicial. In the order constrained statistical inference literature a few orderings that commonly occur are chains, the tree order where $x_{1}\prec x_{2},\ldots,x_{p}$, and the umbrella order where $x_{1}\prec\cdots\prec x_{l-1}\prec x_{l}\succ x_{l+1}\cdots\succ x_{p}$ [49, Sec 2.3]. The umbrella order represents a kind of unimodality where the mode is known to occur at $x_{l}$. Only the umbrella order out of these three orders is not simplicial.

We now pair Lemma 2 with Lemma 1 to find a $\mathcal{V}$-representation of $\mathcal{N}_{<\infty}$. To illustrate, set $\mathcal{P}_{1}=\{1,2\}$ and $\mathcal{P}_{2}=\{1,2,3\}$, both with the standard ordering. The extremal rays of $\mathcal{C}(\mathcal{P}_{1}\times\mathcal{P}_{2})$ are all possible “staircase” matrices, corresponding to upsets of $\mathcal{P}_{1}\times\mathcal{P}_{2}$, that have a staircase of ones in the lower-right corner of the matrix:

[TABLE]

Any connected upset of the poset $\{1,\ldots,m\}$ with the standard ordering has the form $\{a,a+1,\ldots,m\}$ for some $a\in[m]$. Consequently, the extremal rays of $\mathcal{C}(\mathcal{P}_{1})\otimes\mathcal{C}(\mathcal{P}_{2})$ have the form $\big(\sum_{i=a}^{2}{\mathbf{e}}_{i}\big)\big(\sum_{j=b}^{3}{\mathbf{e}}_{i})^{\intercal}$ for $a\in[2]$, $b\in[3]$. These extremal rays are enumerated in the first row of (4) and consist of exactly the “box” matrices that have a rectangle of ones in the lower-right corner. The three matrices in (5) are monotone with respect to the partial order $\mathcal{P}_{1}\times\mathcal{P}_{2}$ but do not possess ND factorizations. In general, the extremal rays of $\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})$ for arbitrary posets $\mathcal{P}_{j}$ have the form $\otimes_{j=1}^{k}(\sum_{i_{j}\in I_{j}}{\mathbf{e}}_{i_{j}})$, where $I_{j}$ is a non-empty, connected upset of $\mathcal{P}_{j}$.

Theorem 3**.**

When $\mathcal{P}=\times_{j=1}^{k}\mathcal{P}_{j}$ the cone $\mathcal{C}(\mathcal{P})$ is equal to $\mathcal{N}_{<\infty}$ if and only if every poset $\mathcal{P}_{j}$ except for one is trivial. A trivial poset has no arrows in its Hasse diagram. If every $\mathcal{P}_{j}=[p_{j}]=\{1,\ldots,p_{j}\}$ is a chain, the cone $\mathcal{N}_{<\infty}$ is generated by $\prod_{j=1}^{k}p_{j}$ extremal rays, while $\mathcal{C}(\mathcal{P})$ has the same number of extremal rays as there are non-empty subsets $\mathcal{A}\subset\mathcal{P}$, where no two elements of $\mathcal{A}$ are comparable. In particular, when $k=2$ the number of extremal rays of $\mathcal{C}(\mathcal{P})$ is ${p_{1}+p_{2}\choose p_{1}}-1$.

Proof.

Let $\mathcal{P}_{j}$ have size $p_{j}$ and have $q_{j}$ connected upsets. The extremal rays of $\mathcal{N}_{<\infty}$ are all extremal rays of $\mathcal{C}(\mathcal{P})$. To see this, observe that every extremal ray of $\mathcal{N}_{<\infty}$ has the form $\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}$ where ${\mathbf{v}}^{(j)}$ is extremal in $\mathcal{C}(\mathcal{P}_{j})$ by Lemma 1. By Lemma 2 each ${\mathbf{v}}^{(j)}$ is associated with a connected upset $\mathcal{U}_{j}\subset\mathcal{P}_{j}$. Any product $\times_{j=1}^{k}\mathcal{U}_{j}$ of connected upsets is a connected upset in $\mathcal{P}$. As $\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}$ is associated with the connected upset $\times_{j=1}^{k}\mathcal{U}_{j}$, another invocation of Lemma 2 shows that $\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}$ is extremal in $\mathcal{C}(\mathcal{P})$.

As a consequence, $\mathcal{N}_{<\infty}=\mathcal{C}(\mathcal{P})$ if and only if these two cones have the same number of extremal rays. Assume that $\mathcal{P}_{1}$ is the only non-trivial poset. Then $\times_{j=1}^{k}\mathcal{P}_{j}$ is a disjoint union of $\prod_{j=2}^{k}p_{j}$ copies of the poset $\mathcal{P}_{1}$. The number of connected upsets of $\mathcal{P}$ is $q_{1}\prod_{j=2}^{k}p_{j}$ and is equal to the number of extremal rays of $\mathcal{C}(\mathcal{P})$. As $p_{j}=q_{j}$ for $j>1$, $q_{1}\prod_{j=2}^{k}p_{j}=\prod_{i=1}^{k}q_{j}$ is equal to the product of the number of connected upsets of each $\mathcal{P}_{j}$, namely the number of extremal rays of $\mathcal{N}_{<\infty}$.

To show the other direction, assume that the inequalities $a_{j}\prec b_{j}$ are present in $\mathcal{P}_{j}$ for $j=1,2$. Letting $x_{j}$ be a maximal element of $\mathcal{P}_{j}$, $j=3,\ldots,k$, the set

[TABLE]

is a connected upset that contains elements of the form $(a_{1},b_{2},x_{3},\ldots,x_{k})$ and $(b_{1},a_{2},x_{3},\ldots,x_{k})$, but does not contain $(a_{1},a_{2},x_{3},\ldots,x_{k})$, implying that it cannot equal a Cartesian product upset $\times_{j=1}^{k}\mathcal{U}_{j}$. The number of connected upsets of $\mathcal{P}$ is at least equal to $\prod_{j=1}^{k}q_{j}+1$ — the number of connected, product upsets plus the non-product upset described above. Therefore $\mathcal{C}(\mathcal{P})$ has more extremal rays than $\prod_{j=1}^{k}q_{j}$, the number of extremal rays of $\mathcal{N}_{<\infty}$.

To prove the last result we count the number of connected upsets of $\times_{j=1}^{k}[p_{j}]$. This is equal to the number of non-empty antichains — sets where no two elements in a set are comparable. There is a bijection between antichains $\mathcal{A}\subset\mathcal{P}$ and connected upsets $\cup_{a\in\mathcal{A}}[a,\infty)$. If $\mathcal{A}=\{a_{1},\ldots,a_{m}\}$ are the minimal elements of $\mathcal{U}$ then $\mathcal{A}$ is an antichain with $\mathcal{U}=\cup_{a\in\mathcal{A}}[a,\infty)$ and conversely any antichain $\mathcal{A}$ produces the upset $\cup_{a\in\mathcal{A}}[a,\infty)$ with minimal elements $\mathcal{A}$. This upset, and in fact any upset, is connected when the poset contains a maximum element, as is the case when $\mathcal{P}=\times_{j=1}^{k}[p_{j}]$ is a product of chains. When $k=2$ the number of antichains is equal to the number of staircase matrices; if $(x_{1},y_{1}),\ldots,(x_{l},y_{l})$ is an antichain then this determines a unique staircase matrix ${\mathbf{M}}$, where $M_{ab}=1$ if and only if there exists an $i=1,\ldots,l$ with $(a,b)\succeq(x_{i},y_{i})$. A staircase matrix is also determined by a path with non-decreasing coordinates in the lattice $\mathbb{N}\times\mathbb{N}$ from $(0,0)$ to $(p_{1},p_{2})$. Such a path must have length $p_{1}+p_{2}$, and out of the $p_{1}+p_{2}$ steps along the path, $p_{1}$ of these steps must be upward in the direction $(1,0)$. The number of such paths is ${p_{1}+p_{2}\choose p_{1}}$. As the path $(0,0)-(1,0)-\cdots-(p_{1},0)-(p_{1},1)-\cdots-(p_{1},p_{2})$ corresponds to the ${\mathbf{0}}$ staircase matrix, or the empty antichain, we must remove this path from consideration. ∎

An implication of Theorem 3 is that if a random tensor ${\mathbf{T}}$ is drawn from a measure that is absolutely continuous with respect to the Lebesgue measure on $\mathcal{C}(\mathcal{P})$ then $\text{Pr}({\mathbf{T}}\in\mathcal{N}_{<\infty})>0$ and $\text{Pr}({\mathbf{T}}\in\mathcal{C}(\mathcal{P})\backslash\mathcal{N}_{<\infty})>0$. The probability that a ${\mathbf{T}}$ drawn from the uniform distribution on the order polytope $\mathcal{O}(\mathcal{P})$ lies in $\mathcal{N}_{<\infty}\cap\mathcal{O}(\mathcal{P})$ when $\mathcal{P}=[m]\times[m]$ is found to equal $50\%,2.38\%$, and $0.004\%$ when $m=2,3,4$ by a volume computation. As the dimension of the matrix grows, the number of extremal rays of $\mathcal{C}(\mathcal{P})$ becomes significantly larger than that of $\mathcal{N}_{<\infty}$, and the probability of observing a random monotone tensor with finite ND rank is small. This suggests that when a data matrix ${\mathbf{T}}={\mathbf{X}}+\bm{\epsilon}$ is observed by taking a finite ND rank matrix ${\mathbf{X}}$ and perturbing it by noise $\bm{\epsilon}$, there is only a small probability that the observed ${\mathbf{T}}$ has finite ND rank. From a parameter estimation perspective this does not pose a serious issue as the observed ${\mathbf{T}}$ can be projected onto $\mathcal{N}_{<\infty}$ or $\mathcal{N}_{\leq r}$. We consider such low-rank approximations in Section 10.

Given a tensor ${\mathbf{T}}$, it must first be ascertained that ${\mathbf{T}}\in\mathcal{N}_{<\infty}$ if ${\mathbf{T}}$ is to have an exact ND factorization. Equations for determining if ${\mathbf{T}}\in\mathcal{N}_{<\infty}$ can be found algorithmically by converting the $\mathcal{V}$-representation of $\mathcal{N}_{<\infty}$ in Lemma 1 into a $\mathcal{H}$-representation. In the next section we provide simple equations for testing membership when every $\mathcal{P}_{i}$ except for one has a tree structure. Below, the $24$ equations of $\mathcal{N}_{<\infty}\subset\mathbb{R}_{+}^{3\times 3}$ are computed in Polymake [27] for the case where $\mathcal{P}_{i}=\{a^{(i)},b^{(i)},c^{(i)}\}$, $i=1,2$ is given by the collider structure specified in Definition 4, and the third row and column respectively correspond to the maximum elements $c^{(1)},c^{(2)}$.

[TABLE]

The first three rows of equations correspond to simple positivity and monotonicity constraints. The next two rows of equations have four non-zero variables and will appear in the next section as the type of equations required to cut-out $\mathcal{N}_{<\infty}$ when the posets are trees. The final four rows of equations are more exotic as they can be written as $\text{tr}({\mathbf{C}}^{\intercal}{\mathbf{T}})\leq 0$ where the matrix ${\mathbf{C}}$ has a rank greater than one. For example, the last two inequalities have the corresponding ${\mathbf{C}}$ matrices

[TABLE]

These matrices do not lie in $\otimes_{j=1}^{2}\mathcal{C}(\mathcal{P}_{j})^{*}$, for if they did they would be extremal in $\otimes_{j=1}^{2}\mathcal{C}(\mathcal{P}_{j})^{*}$ as they are normal vectors of facets of $\otimes_{j=1}^{2}\mathcal{C}(\mathcal{P}_{j})$, and by Lemma 1 they would be rank one matrices.

6 Reduction of the ND rank to Nonnegative Rank

For invertible linear maps ${\mathbf{A}}^{(j)}$, we define the tensor product of these linear maps as the linear map ${\mathbf{A}}^{(1)}\otimes\cdots\otimes{\mathbf{A}}^{(k)}:\mathbb{R}^{p_{1}\times\cdots\times p_{k}}\rightarrow\mathbb{R}^{p_{1}\times\cdots\times p_{k}}$ that is given in coordinates by

[TABLE]

On rank-one tensors this linear map is given by the simple formula $(\otimes_{j=1}^{k}{\mathbf{A}}^{(j)})(\otimes_{j=1}^{k}{\mathbf{v}}^{(j)})=\otimes_{j=1}^{k}{\mathbf{A}}^{(j)}{\mathbf{v}}^{(j)}$. The idea in this section is to choose the maps ${\mathbf{A}}^{(j)}$ so that the order cone $\mathcal{C}(\mathcal{P}_{j})$ is mapped bijectively onto the positive orthant $\mathbb{R}^{p_{j}}_{+}$ by ${\mathbf{A}}^{(j)}$. As the tensor product map preserves rank-one tensors, the problem of determining the ND rank can be transformed into a problem of determining the nonnegative rank whenever such ${\mathbf{A}}^{(j)}$ exist.

Theorem 4**.**

Assume that $\mathcal{C}(\mathcal{P}_{j})$ is simplicial for every $j=1,\ldots,k$ and let ${\mathbf{A}}^{(j)}\in\mathrm{GL}(p_{j})$ be an invertible linear map where ${\mathbf{A}}^{(j)}(\mathcal{C}(\mathcal{P}_{j}))=\mathbb{R}^{p_{j}}_{+}$. If ${\mathbf{T}}\in\mathcal{N}_{<\infty}$ then $\mathrm{NDrank}({\mathbf{T}})=\mathrm{rank}_{+}\big(\otimes_{j=1}^{k}{\mathbf{A}}^{(j)}({\mathbf{T}})\big)$.

Proof.

As $\mathcal{C}(\mathcal{P}_{j})$ is simplicial it is the conical hull of ${\mathbf{v}}^{(1)},\ldots,{\mathbf{v}}^{(p_{j})}$. We can define ${\mathbf{A}}^{(j)}$ to be the linear map that sends ${\mathbf{v}}^{(i)}$ to ${\mathbf{e}}_{i}$ and thus maps $\mathcal{C}(\mathcal{P}_{j})$ onto $\mathbb{R}^{p_{j}}_{+}$. If ${\mathbf{T}}=\sum_{i=1}^{r}\otimes_{j=1}^{k}{\mathbf{v}}^{(ij)}$ is a rank-$r$ ND factorization of ${\mathbf{T}}$ then $\otimes_{j=1}^{k}{\mathbf{A}}^{(j)}({\mathbf{T}})=\sum_{i=1}^{r}\otimes_{j=1}^{k}{\mathbf{A}}^{(j)}{\mathbf{v}}^{(ij)}$ is a nonnegative rank-$r$ tensor factorization of $\otimes_{j=1}^{k}{\mathbf{A}}^{(j)}({\mathbf{T}})$. Conversely, if $\otimes_{j=1}^{k}{\mathbf{A}}^{(j)}({\mathbf{T}})=\sum_{i=1}^{r}\otimes_{j=1}^{k}{\mathbf{w}}^{(ij)}$ is a rank-$r$ nonnegative factorization of $\otimes_{j=1}^{k}{\mathbf{A}}^{(j)}({\mathbf{T}})$ then $\sum_{i=1}^{r}\otimes_{j=1}^{k}({\mathbf{A}}^{(j)})^{-1}{\mathbf{w}}^{(ij)}$ is a rank-$r$ ND factorization of ${\mathbf{T}}$. ∎

By Lemma 2, $\mathcal{C}(\mathcal{P}_{j})$ is simplicial when the Hasse diagram of $\mathcal{P}_{j}$ is a union of trees. One special case where the above theorem applies is when $\mathcal{P}=\times_{j=1}^{k}[p_{j}]$ under the standard ordering; the problem of finding an exact ND factorization can be converted into the problem of finding a nonnegative factorization. The matrices ${\mathbf{A}}^{(j)}$ required for this conversion are related to Möbius inversions. The Möbius tranform ${\mathbf{M}}_{\mathcal{P}}:\mathbb{R}^{\mathcal{P}}\rightarrow\mathbb{R}^{\mathcal{P}}$ of a poset $\mathcal{P}$ is defined as the linear map with

[TABLE]

where $\{{\mathbf{e}}_{x}\}_{x\in\mathcal{P}}$ is the standard basis of $\mathbb{R}^{\mathcal{P}}$.

Lemma 3**.**

If $\mathcal{P}$ has connected components $\mathcal{S}_{1},\ldots,\mathcal{S}_{m}$ then ${\mathbf{M}}_{\mathcal{P}}=\oplus_{i=1}^{m}{\mathbf{M}}_{\mathcal{S}_{i}}$. Moreover, if $\mathcal{P}=\times_{j=1}^{k}\mathcal{P}_{j}$ then ${\mathbf{M}}_{\mathcal{P}}=\otimes_{j=1}^{k}{\mathbf{M}}_{\mathcal{P}_{j}}$.

Proof.

We first prove that ${\mathbf{M}}_{\mathcal{P}}=\oplus_{i=1}^{m}{\mathbf{M}}_{\mathcal{S}_{i}}$, meaning that ${\mathbf{M}}_{\mathcal{P}}$ is a block-diagonal matrix with $m$ diagonal blocks that equal the ${\mathbf{M}}_{\mathcal{S}_{i}}$s. If $x\in\mathcal{S}_{i}$ then ${\mathbf{M}}_{\mathcal{P}}{\mathbf{e}}_{x}=\sum_{y\in\mathcal{P}:x\preceq y}{\mathbf{e}}_{y}=\sum_{y\in\mathcal{S}_{i}:x\preceq y}{\mathbf{e}}_{y}={\mathbf{M}}_{\mathcal{S}_{i}}{\mathbf{e}}_{x}$, where we have used the fact that the only elements in $\mathcal{P}$ that are comparable to $x$ are in $\mathcal{S}_{i}$. This proves the first statement. The second statement follows from

[TABLE]

where ${\mathbf{e}}_{(x_{1},\ldots,x_{k})}=\otimes_{j=1}^{k}{\mathbf{e}}_{x_{j}}$ is a rank-one, standard basis element of $\mathbb{R}^{\mathcal{P}}=\otimes_{j=1}^{k}\mathbb{R}^{\mathcal{P}_{j}}$. That ${\mathbf{M}}_{\mathcal{P}}=\otimes_{j=1}^{k}{\mathbf{M}}_{\mathcal{P}_{j}}$ is a consequence of the tensor product definition for linear maps provided in (6). ∎

Lemma 4**.**

If $\mathcal{P}$ is simplicial, ${\mathbf{M}}_{\mathcal{P}}^{-1}$ bijectively maps $\mathcal{C}(\mathcal{P})$ onto $\mathbb{R}^{p}_{+}$.

Proof.

When $\mathcal{C}(\mathcal{P})$ is simplicial every upset $\mathcal{U}$ of $\mathcal{P}$ has the form $\mathcal{U}=[x,\infty)$ by the proof in Theorem 2. The column corresponding to ${\mathbf{e}}_{x}$ in ${\mathbf{M}}_{\mathcal{P}}$ is equal to $\sum_{y\in[x,\infty)}{\mathbf{e}}_{y}$, which is one of the extremal rays of $\mathcal{C}(\mathcal{P})$ by Lemma 2. Thus, the columns of ${\mathbf{M}}_{\mathcal{P}}$ are exactly the extremal rays of $\mathcal{C}(\mathcal{P})$ and so ${\mathbf{M}}_{\mathcal{P}}$ maps $\mathbb{R}^{p}_{+}$ bijectively onto $\mathcal{C}(\mathcal{P})=\text{cone}([{\mathbf{M}}_{\mathcal{P}}]_{\cdot 1},\ldots,[{\mathbf{M}}_{\mathcal{P}}]_{\cdot p})$. ∎

Lemma 5**.**

If $\mathcal{C}(\mathcal{P})$ is simplicial, the matrix ${\mathbf{M}}_{\mathcal{P}}^{-1}$ has entries

[TABLE]

When $\mathcal{P}=[p]=\{1,\ldots,p\}$ is a chain the above formula simplifies to a Toeplitz matrix of the form:

[TABLE]

where the columns and rows are ordered with respect to the basis ${\mathbf{e}}_{1},\ldots,{\mathbf{e}}_{p}$.

Proof.

We compute

[TABLE]

The last equality follows from the assumed tree structure of $\mathcal{P}$, where the upset $[y,\infty)$ consists of $y$ along with the upsets of all elements $z$ that cover $y$. If $z_{1},z_{2}$ both cover $y$ then $[z_{1},\infty)\cap[z_{2},\infty)=\emptyset$ since $\mathcal{P}$ has no colliders. This ensures that there are no repeats of basis elements in the double sum above. ∎

When each $\mathcal{P}_{j}=[p_{j}]=\{1,\ldots,p_{j}\}$ is a chain, Lemma 4 implies that ${\mathbf{M}}_{\mathcal{P}_{j}}^{-1}$ bijectively maps $\mathcal{C}(\mathcal{P}_{j})$ onto $\mathbb{R}^{p_{j}}_{+}$, while Lemma 3 further implies that $\otimes_{j=1}^{k}{\mathbf{M}}_{\mathcal{P}_{j}}^{-1}$ maps $\mathcal{N}_{<\infty}$ bijectively onto $\otimes_{j=1}^{k}\mathbb{R}^{p_{j}}_{+}$. An explicit formula for ${\mathbf{M}}_{\mathcal{P}_{j}}^{-1}$ is provided in (7). This provides a certificate of whether ${\mathbf{T}}$ is in $\mathcal{N}_{<\infty}$, since this occurs if and only if the entries of $(\otimes_{j=1}^{k}{\mathbf{M}}_{\mathcal{P}_{j}}^{-1})({\mathbf{T}})$ are nonnegative. The $(i_{1},\ldots,i_{k})$ entry of this tensor is

[TABLE]

where we use the convention that $T_{i_{1}-j_{1},\ldots,i_{k}-j_{k}}=0$ whenever $i_{l}-j_{l}=0$ for any $l$. Viewing ${\mathbf{T}}$ as a function, we can define the differencing operator along the $j$th mode as $(\Delta^{(j)}{\mathbf{T}})(i_{1},\ldots,i_{k})=T(i_{1},\ldots,i_{j},\ldots,i_{k})-T(i_{1},\ldots,i_{j}-1,\ldots,i_{k})$. From (8) the tensor $(\otimes_{j=1}^{k}{\mathbf{M}}_{\mathcal{P}_{j}}^{-1})({\mathbf{T}})$ can be seen to be equal to the tensor $\Delta^{(1)}\cdots\Delta^{(k)}{\mathbf{T}}$. In the matrix setting, with $p_{1}=2$, $p_{2}=3$ we obtain

[TABLE]

Equations of the above form were highlighted in [45], where monotonicity for functions defined on a Cartesian product of intervals were examined. The notion of a ${\mathbf{1}}_{k}$ increasing tensor in [45] is equivalent to the present notion of a finite ND rank tensor. Notice that similar expressions appeared in the previous section for the $\mathcal{H}$-representation for a product of collider posets.

We now discuss a probabilistic interpretation of the Möbius transform of the poset $\mathcal{P}=\times_{j=1}^{k}[p_{j}]$ and Theorem 4. Suppose that ${\mathbf{R}}$ is a nonnegative tensor with entries summing to one. The entry $R_{i_{1}\ldots i_{k}}$ can be interpreted as a probability, where ${\mathbf{R}}$ is a probability mass function (PMF) on $\mathcal{P}=\times_{j=1}^{k}[p_{j}]$. Any rank-$r$ nonnegative tensor factorization of ${\mathbf{R}}$ can be written in the form $\sum_{i=1}^{r}\lambda_{i}\otimes_{j=1}^{k}{\mathbf{q}}^{(ij)}$ with vectors $\bm{\lambda}\in\Delta_{r-1}$ and ${\mathbf{q}}^{(ij)}\in\Delta_{p_{j}-1}$ all residing in probability simplices. This means that ${\mathbf{R}}$ is a mixture of $r$ probability mass functions $\otimes_{j=1}^{k}{\mathbf{q}}^{(ij)}$ on $\mathcal{P}$, with respective mixture weights $\lambda_{i}$ [22, Ch 4]. The tensor ${\mathbf{M}}_{\mathcal{P}}({\mathbf{R}})$ is exactly the multivariate conditional distribution function (CDF) of ${\mathbf{R}}$, meaning that if ${\mathbf{X}}\in\mathcal{P}$ is a random variable with distribution given by ${\mathbf{R}}$ then ${\mathbf{M}}_{\mathcal{P}}({\mathbf{R}})_{i_{1}\ldots i_{k}}=\text{Pr}({\mathbf{X}}\preceq(i_{1},\ldots,i_{k}))$. If ${\mathbf{F}}_{i}\in\mathcal{N}_{1}$ is the CDF corresponding to $\otimes_{j=1}^{k}{\mathbf{q}}^{(ij)}$ then ${\mathbf{M}}_{\mathcal{P}}({\mathbf{R}})=\sum_{i=1}^{r}\lambda_{i}{\mathbf{F}}_{i}$ has an ND rank of $r$. In summary, a probability distribution has a PMF with a nonnegative rank of $r$ if and only if it has a CDF with a ND rank of $r$.

To conclude this section we provide the $\mathcal{H}$-representation for $\mathcal{N}_{<\infty}$ when all but one cone $\mathcal{C}(\mathcal{P}_{j})$ is simplicial. The following theorem generalizes equation (8) beyond the case where each $\mathcal{P}_{j}$ is a chain. Furthermore, the left-hand side of equation (9) is the formula for obtaining a PMF from the CDF ${\mathbf{T}}$ when each $\mathcal{P}_{j}$ is a chain.

Theorem 5**.**

Assume that $\mathcal{P}_{j}$ has no colliders for all but one $j=1,\ldots,k$. Define the poset $\mathcal{P}_{j}^{\prime}=\{0_{j}\}\cup\mathcal{P}_{j}$ where $0_{j}\prec x$ for all $x\in\mathcal{P}_{j}$. For every pair of elements $(x_{j},y_{j})$ in $\mathcal{P}_{j}$ with $y_{j}$ covering $x_{j}$ define the inequality

[TABLE]

where $\delta(i_{j})=0$ if $i_{j}=y_{j}$ and $1$ when $i_{j}=x_{j}$. The set of finite ND rank tensors with respect to $\times_{j=1}^{k}\mathcal{P}_{j}$ is the cone defined by the intersection of all of the halfspaces (9). Each hyperplane in (9) supports a facet of $\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})$.

Proof.

By Theorem 1 $\big(\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})\big)^{*}=\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})^{*}$. The facets of $\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})$ are supported by hyperplanes corresponding to the extremal rays of $\big(\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})\big)^{*}$. Lemma 1 shows that the extremal rays of $\big(\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})\big)^{*}$ equal $\otimes_{j=1}^{k}{\mathbf{h}}^{(j)}$ where each ${\mathbf{h}}^{(j)}$ is an extremal ray of $\mathcal{C}(\mathcal{P}_{j})^{*}$. The set $\mathcal{C}(\mathcal{P}_{j})$ is the intersection of the halfspaces of the form $\tilde{H}_{xy}=\{{\mathbf{f}}:f_{x}\leq f_{y}\}$ for $x\preceq y$ and $x,y\in\mathcal{P}_{j}^{\prime}$. In the case where $x=0$ we define $\tilde{H}_{0y}=\{{\mathbf{f}}:0\leq f_{y}\}$ and we let $H_{xy}$ be the hyperplane associated with the halfspace $\tilde{H}_{xy}$. The intersection $H_{xy}\cap\mathcal{C}(\mathcal{P}_{j})$ is only $(p_{j}-1)$-dimensional when $x\lessdot y$ in $\mathcal{P}_{j}^{\prime}$, as if $x\prec z\prec y$ then $f_{x}=f_{z}=f_{y}$ for all ${\mathbf{f}}$ in $H_{xy}$ (where we have also defined $f_{0}\coloneq 0$) and $H_{xy}\cap\mathcal{C}(\mathcal{P}_{j})$ has dimension less than $p_{j}-1$. Conversely, if $x\lessdot y$ then $H_{ab}\cap\mathcal{C}(\mathcal{P}_{j})$ has dimension $p_{j}-1$.

To complete the proof, let ${\mathbf{h}}^{(j)}={\mathbf{e}}_{y_{j}}^{*}-{\mathbf{e}}_{x_{j}}^{*}$ for $x_{j}\lessdot y_{j}$ in $\mathcal{P}_{j}$, where ${\mathbf{e}}_{0}\coloneq{\mathbf{0}}$ and each ${\mathbf{e}}_{x}^{*}$ is the dual vector for ${\mathbf{e}}_{x}$ with respect to the standard bases $\{{\mathbf{e}}_{x}\}_{x\in\mathcal{P}_{j}}$ of each $\mathcal{P}_{j}$. The hyperplane defining the facet $\otimes_{j=1}^{k}{\mathbf{h}}^{(j)}$ is

[TABLE]

∎

From this theorem we see that the matrix appearing in Table 1 in the selenium example does not possess an exact ND decomposition. This is because the constraints, $t_{12}\leq t_{32}$, $t_{11}-t_{12}-t_{31}+t_{32}\geq 0$, and $t_{13}-t_{14}-t_{33}+t_{34}\geq 0$, all of which are necessary for a ND factorization to exist, do not hold.

7 Maximum and Typical ND Ranks

7.1 The Maximum ND Rank

The maximum (finite) ND rank associated with the order cones $\mathcal{C}(\mathcal{P}_{j})$, $j=1,\ldots,k$ is defined as

[TABLE]

In this section the maximum ND rank is found for certain order cones in the matrix setting. Even in the matrix case, finding the maximum ND rank for general order cones is a challenging problem. In comparison, the maximum nonnegative rank, is known to equal $\prod_{j=1}^{k-1}p_{j}$ for tensors in $\otimes_{j=1}^{k}\mathbb{R}^{p_{j}}_{+}$ where $p_{1}\leq\cdots\leq p_{k}$ [54]. An analogous upper bound on the maximum ND rank is found below.

Lemma 6**.**

Let $\mathcal{C}(\mathcal{P}_{j})$ have $q_{j}$ extremal rays where $q_{1}\leq\cdots\leq q_{k}$. The maximum ND rank in $\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})$ is at most $\prod_{j=1}^{k-1}q_{j}$.

Proof.

Let ${\mathbf{T}}\in\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})$ have the ND factorization $\sum_{i=1}^{r}\otimes_{j=1}^{k}{\mathbf{v}}^{(ij)}$. Letting ${\mathbf{w}}^{(1j)},\ldots,{\mathbf{w}}^{(q_{j}j)}$ be the extremal rays of $\mathcal{C}(\mathcal{P}_{j})$, there exists $\alpha_{ijl}$ with ${\mathbf{v}}^{(ij)}=\sum_{l=1}^{q_{j}}\alpha_{ijl}{\mathbf{w}}^{(lj)}$. We have that

[TABLE]

which is a sum of $\prod_{j=1}^{k-1}q_{j}$ rank one terms in $\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})$. ∎

If there are a large number of extremal rays in the given order cones the following lemma can provide a tighter upper bound on the maximum ND rank.

Lemma 7**.**

The maximum ND rank is at most $\prod_{j=1}^{k}p_{j}$, where $\mathcal{C}(\mathcal{P}_{j})\subset\mathbb{R}^{p_{j}}$.

Proof.

By Caratheodory’s theorem, for any tensor ${\mathbf{T}}\in\otimes_{j=1}^{k}\mathcal{C}(\mathcal{P}_{j})=\text{cone}(\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}:{\mathbf{v}}^{(j)}\in\mathcal{C}(\mathcal{P}_{j}))$ the number of elements in $\{\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}:{\mathbf{v}}^{(j)}\in\mathcal{C}(\mathcal{P}_{j})\}$ needed to write ${\mathbf{T}}$ as a conic combination is at most $\prod_{j=1}^{k}p_{j}$. ∎

To find a lower bound on the maximum ND rank we use an analogue to the nested cone condition [29, Sec 2.1.1] for nonnegative matrix factorizations.

Lemma 8**.**

The ND rank of a matrix ${\mathbf{T}}\in\mathcal{N}_{<\infty}$ is at least equal to the minimum number of elements ${\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{s}\in\mathcal{C}(\mathcal{P}_{1})$ needed for ${\mathbf{T}}(\mathcal{C}(\mathcal{P}_{2})^{*})\subseteq\mathrm{cone}({\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{s})$.

Proof.

Suppose that ${\mathbf{T}}=\sum_{i=1}^{r}{\mathbf{a}}_{i}\otimes{\mathbf{b}}_{i}$ is a ND factorization. Then ${\mathbf{T}}(\mathcal{C}(\mathcal{P}_{2})^{*})\subseteq\mathrm{cone}({\mathbf{a}}_{1},\ldots,{\mathbf{a}}_{r})$ since ${\mathbf{T}}(\bm{\beta})=\sum_{i=1}^{r}\bm{\beta}({\mathbf{b}}_{i}){\mathbf{a}}_{i}\in\mathrm{cone}({\mathbf{a}}_{1},\ldots,{\mathbf{a}}_{r})$ whenever $\bm{\beta}\in\mathcal{C}(\mathcal{P}_{2})^{*}$. The number $r$ is at least equal to $s$ by assumption. ∎

The upper bound in Lemma 6 is exact when one of the constituent order cones is simplicial. The posets occurring in many applications will satisfy this condition; the selenium example in Section 4 satisfies this condition and has a maximum rank of four.

Theorem 6**.**

If $\mathcal{C}(\mathcal{P}_{j})$ has $q_{j}$ extremal rays, $j=1,2$, and $\mathcal{C}(\mathcal{P}_{2})$ is simplicial then the maximum ND rank of a matrix in $\mathcal{N}_{<\infty}$ is $\min(q_{1},q_{2})$.

Proof.

Without loss of generality it may be assumed that $\mathcal{C}(\mathcal{P}_{2})=\mathbb{R}^{q_{2}}_{+}$ by the affine-invariance of the ND rank (Theorem 4). As $\mathcal{C}(\mathcal{P}_{2})$ is simplicial $p_{2}=q_{2}$. Let ${\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{q_{1}}$ be the extremal rays of $\mathcal{C}(\mathcal{P}_{1})$. If $q_{2}\leq p_{1}\leq q_{1}$ then any full rank matrix ${\mathbf{T}}$ will have $q_{2}=\text{rank}({\mathbf{T}})\leq\text{NDrank}({\mathbf{T}})\leq q_{2}$ by Lemma 6. If $p_{1}\leq q_{2}\leq q_{1}$ take ${\mathbf{T}}=\sum_{i=1}^{q_{2}}{\mathbf{v}}_{i}\otimes{\mathbf{e}}_{i}$. Then ${\mathbf{T}}(\mathcal{C}(\mathcal{P}_{2})^{*})={\mathbf{T}}(\mathbb{R}^{q_{2}}_{+})=\text{cone}({\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{q_{2}})$. If $\text{cone}({\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{q_{2}})\subseteq\text{cone}({\mathbf{a}}_{1},\ldots,{\mathbf{a}}_{t})$ for some ${\mathbf{a}}_{i}\in\mathcal{C}(\mathcal{P}_{1})$ then vectors proportional to ${\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{q_{2}}$ must appear in the various ${\mathbf{a}}_{i}$s as the ${\mathbf{v}}_{i}$s are extremal. Using Lemma 8 the ND rank of ${\mathbf{T}}$ is $q_{2}$. Finally, when $p_{1}\leq q_{1}\leq q_{2}$, taking ${\mathbf{T}}=\sum_{i=1}^{q_{1}}{\mathbf{v}}_{i}\otimes{\mathbf{e}}_{i}$, a similar argument using ${\mathbf{T}}(\mathcal{C}(\mathcal{P}_{2})^{*})={\mathbf{T}}(\mathbb{R}^{q_{2}}_{+})=\text{cone}({\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{q_{1}})=\mathcal{C}(\mathcal{P}_{1})$ shows that the ND rank of ${\mathbf{T}}$ is $q_{1}$. ∎

Of note in the above theorem is that the maximum ND rank of a matrix in $\mathbb{R}^{p_{1}\times p_{2}}$ can potentially be much larger than the usual maximum matrix rank of $\min(p_{1},p_{2})$.

The prototypical example of a poset with a non-simplicial order cone is a collider (Definition 4). The next theorem shows that when the column and row posets of a matrix both have the form of a collider the maximum ND rank can be larger than even $\max(p_{1},p_{2})$.

Theorem 7**.**

Let $\mathcal{P}=\{x_{1},\ldots,x_{p}\}$ where $x_{i}\prec x_{p}$ for all $i<p$. The maximum rank in $\mathcal{C}(\mathcal{P})\otimes\mathcal{C}(\mathcal{P})\subset\mathbb{R}^{p\times p}$ is at least $2(p-1)$. When $p=3$ the maximum rank is equal to $4$.

Proof.

Let $p\geq 3$ so that $\mathcal{P}$ contains a collider. The extremal rays ${\mathbf{v}}_{1},\ldots,{\mathbf{v}}_{2^{p-1}}$, in $\mathcal{C}(\mathcal{P})$ can be taken to have the form $({\mathbf{w}},1)$ where ${\mathbf{w}}\in\mathbb{R}^{p-1}$ is any one of the $2^{p-1}$ vectors that have entries equal to either zero or one. Hence, $\mathcal{C}(\mathcal{P})$ is equal to the homogenization [59, Sec 1.5] of a hypercube. The facets of $\mathcal{C}(\mathcal{P})$ each contain $2^{p-2}$ extremal rays and there are $2(p-1)$ facets; one facet for each equality $0=f_{x_{i}}$, and one facet for $f_{x_{i}}=f_{x_{p}}$, where $i<p$. The facet $\{{\mathbf{f}}:f_{x_{i}}=0\}\cap\mathcal{C}(\mathcal{P})$ contains the extremal rays of the form $({\mathbf{w}},1)$ with ${\mathbf{w}}_{i}=0$, while the facet $\{{\mathbf{f}}:f_{x_{i}}=f_{x_{p}}\}\cap\mathcal{C}(\mathcal{P})$ contains the remaining extremal rays $({\mathbf{w}},1)$ with $w_{i}=1$. These two facets are opposite to each other and partition the set of extremal rays.

The goal of the remainder of the proof to construct a matrix ${\mathbf{T}}$ that has an ND rank bounded below by $2(p-1)$. Choose a facet $F$ of $\mathcal{C}(\mathcal{P})$ that has $F^{\prime}$ as its opposite facet. Let $J$ be an index set of size less than $p-1$ where ${\mathbf{v}}_{j}\in F$ for every $j\in J$. The set $\text{cone}({\mathbf{v}}_{j}:j\in J)$ has a dimension that is less than $p-1$. If we remove all such cones from $F$ we are left with a non-empty set

[TABLE]

because $F$ has dimension $p-1$. Choose a point $\sum_{i\in I_{F}}\alpha_{i}{\mathbf{v}}_{i}$ in the above set, where the index subset has $|I_{F}|=p-1$ and ${\mathbf{v}}_{i}\in F$ for every $i\in I_{F}$. Repeat the above process for the facet $F^{\prime}$ to obtain an analogous point $\sum_{i\in I_{F^{\prime}}}\alpha_{i}{\mathbf{v}}_{i}$. We define the matrix ${\mathbf{T}}=\sum_{i\in I}\alpha_{i}{\mathbf{v}}_{i}\otimes{\mathbf{v}}_{i}$ where $I=I_{F}\cup I_{F^{\prime}}$.

Assume that ${\mathbf{T}}=\sum_{i=1}^{r}{\mathbf{a}}_{i}\otimes{\mathbf{b}}_{i}$ is a minimum ND rank decomposition, with ${\mathbf{a}}_{i}=\sum_{j=1}^{2^{p-1}}\lambda_{ij}{\mathbf{v}}_{j}$ and ${\mathbf{b}}_{i}=\sum_{l=1}^{2^{p-1}}\mu_{il}{\mathbf{v}}_{l}$, $\lambda_{ij},\mu_{il}\geq 0$. For any pair of extremal rays ${\mathbf{v}}_{l}\neq{\mathbf{v}}_{j}$ choose a facet $G$ with corresponding dual vector ${\mathbf{h}}$ so that ${\mathbf{v}}_{l}\in G$ but ${\mathbf{v}}_{j}\notin G$. Let ${\mathbf{h}}^{\prime}$ be the dual vector corresponding to the facet $G^{\prime}$ opposite of $G$. We compute

[TABLE]

The first equality follows because every extremal ray in $\mathcal{C}(\mathcal{P})$ must be in only one of the facets $G$ or $G^{\prime}$. The inequality is a result of the middle expression being a sum of nonnegative terms. By construction $({\mathbf{h}}^{\intercal}{\mathbf{v}}_{j})({\mathbf{v}}_{l}^{\intercal}{\mathbf{h}}^{\prime})>0$, which implies that $\lambda_{ij}\mu_{il}=0$ whenever $j\neq l$. We conclude that ${\mathbf{a}}_{i}\otimes{\mathbf{b}}_{i}\,\propto\,{\mathbf{v}}_{j_{i}}\otimes{\mathbf{v}}_{j_{i}}$ for some index $j_{i}$, with ${\mathbf{T}}=\sum_{i=1}^{r}c_{i}{\mathbf{v}}_{j_{i}}\otimes{\mathbf{v}}_{j_{i}}$ and $c_{i}>0$. Let ${\mathbf{h}}$ be a dual vector corresponding to the hyperplane that supports the facet $F^{\prime}$. Then

[TABLE]

where without loss of generality it is assumed that ${\mathbf{h}}$ is scaled so that ${\mathbf{h}}^{\intercal}{\mathbf{v}}_{i}=1$ for ${\mathbf{v}}_{i}\in F$. By how the vector $\sum_{i\in I_{F}}\alpha_{i}{\mathbf{v}}_{i}$ was constructed, there must exist at least $p-1$ terms in the sum $\sum_{i:{\mathbf{v}}_{j_{i}}\in F}c_{i}{\mathbf{v}}_{j_{i}}$. A symmetric argument shows that the set of indices $\{i:{\mathbf{v}}_{j_{i}}\in F^{\prime}\}$ has a size of at least $p-1$. As the facets $F$ and $F^{\prime}$ do not share any extremal rays, we conclude that the sum $\sum_{i=1}^{r}c_{i}{\mathbf{v}}_{j_{i}}\otimes{\mathbf{v}}_{j_{i}}$ must include at least $2(p-1)$ terms and so $r\geq 2(p-1)$.

In the case when $p=3$, Lemma 6 shows that the maximum rank is also at most $4=2(3-1)$, which is the number of vertices in a two-dimensional cube. ∎

The problem of finding the maximum rank of tensor decompositions for arbitrary order cones and for more general polyhedral cones remains an open question.

7.2 A Matrix Tri-Factorization Formulation

A succinct representation of a rank-$r$, nonnegative matrix factorization of ${\mathbf{T}}\in\mathbb{R}^{p_{1}\times p_{2}}_{+}$ is ${\mathbf{T}}={\mathbf{A}}^{(1)}({\mathbf{A}}^{(2)})^{\intercal}$ where ${\mathbf{A}}^{(j)}\in\mathbb{R}_{+}^{p_{j}\times r}$. If ${\mathbf{T}}=\sum_{i=1}^{r}{\mathbf{a}}^{(i1)}\otimes{\mathbf{a}}_{2}^{(i2)}$ is an ND factorization we can likewise write this as ${\mathbf{T}}={\mathbf{A}}^{(1)}({\mathbf{A}}^{(2)})^{\intercal}$ where ${\mathbf{A}}^{(j)}\in\mathbb{R}^{p_{j}\times r}$ has columns that are given by the ${\mathbf{a}}^{(ij)}$s for $i=1,\ldots,r$. Let ${\mathbf{V}}^{(j)}\in\mathbb{R}^{p_{i}\times q_{j}}$ be a matrix that has columns that are equal to the $q_{j}$ extremal rays in $\mathcal{C}(\mathcal{P}_{j})$. As each ${\mathbf{a}}^{(ij)}$ is in the conical hull of the columns of ${\mathbf{V}}^{(j)}$ there exist nonnegative matrices ${\mathbf{H}}^{(i)}\in\mathbb{R}_{+}^{q_{i}\times r}$ where ${\mathbf{A}}^{(j)}={\mathbf{V}}^{(j)}{\mathbf{H}}^{(j)}$ for $j=1,2$. Defining ${\mathbf{H}}={\mathbf{H}}^{(1)}({\mathbf{H}}^{(2)})^{\intercal}\in\mathbb{R}^{q_{1}\times q_{2}}_{+}$ the ND factorization can be written as

[TABLE]

This leads to an observation connecting nondecreasing and nonnegative ranks:

Lemma 9**.**

The nondecreasing rank of a matrix ${\mathbf{T}}\in\mathcal{N}_{<\infty}$ is the smallest nonnegative rank of a matrix ${\mathbf{H}}\in\mathbb{R}^{q_{1}\times q_{2}}$ that satisfies the equation ${\mathbf{T}}={\mathbf{V}}^{(1)}{\mathbf{H}}({\mathbf{V}}^{(2)})^{\intercal}$, where ${\mathbf{V}}^{(j)}\in\mathbb{R}^{p_{j}\times q_{j}}$ are fixed matrices with columns that are equal to the $q_{j}$ extremal rays of $\mathcal{C}(\mathcal{P}_{j})$.

The above tri-factorization formulation is related to nonnegative matrix tri-factorizations [56], as well as convex NMF [21], which specifies a fixed dictionary of vectors to be used in the factorization.

As an example of this factorization, the matrix attaining the maximum rank appearing in the proof of Theorem 7 for $p=3$ has the following representation

[TABLE]

where ${\mathbf{H}}$ has a nonnegative rank of four and the rows of ${\mathbf{V}}^{(1)}={\mathbf{V}}^{(2)}$ are equal to the four extremal rays of the collider $\mathcal{C}(\mathcal{P}_{j})$ outlined in Lemma 2.

7.3 Typical ND Ranks

A typical nondecreasing rank is defined as a number $r$ such that $\mathcal{N}_{r}$ has a non-empty interior. The probabilistic interpretation of an typical ND rank is that if a matrix ${\mathbf{X}}$ is drawn from a distribution supported on $\mathcal{N}_{<\infty}$ that has a density with respect to the Lebesgue measure then $\text{Pr}(\text{NDrank}({\mathbf{X}})=r)>0$ when $r$ is a typical rank. In this section we determine the typical ranks for the two matrix settings outlined in Theorems 6 and 7, as well as for any other setting where the maximum ND rank is known.

Theorem 8**.**

The maximum ND rank is always a typical ND rank.

Proof.

Let $m$ be the maximum ND rank and choose a ${\mathbf{T}}$ with $\text{NDrank}({\mathbf{T}})=m$. Let $B_{n^{-1}}({\mathbf{T}})$ be an open ball of radius $n^{-1}$ centered at ${\mathbf{T}}$. The set $B_{n^{-1}}({\mathbf{T}})\cap\mathcal{N}_{<\infty}$ has a non-empty interior. If ${\mathbf{T}}\in\text{int}(\mathcal{N}_{<\infty})$ this is immediate. Otherwise, if ${\mathbf{T}}$ is in the boundary of $\mathcal{N}_{<\infty}$ take an ${\mathbf{R}}\in\text{int}(\mathcal{N}_{<\infty})$, which exists because $\mathcal{N}_{<\infty}$ is proper, and consider the line segment $\lambda{\mathbf{T}}+(1-\lambda){\mathbf{R}}$, $\lambda\in[0,1]$. For all $\lambda\in[0,1)$ the point $\lambda{\mathbf{T}}+(1-\lambda){\mathbf{R}}$ is in $\text{int}(\mathcal{N}_{<\infty})$, and for large enough $\lambda$ this point is also in $B_{n^{-1}}({\mathbf{T}})$. this provides a point in $B_{n^{-1}}({\mathbf{T}})\cap\text{int}(\mathcal{N}_{<\infty})$. Now if the maximum ND rank was not a typical rank, the set $B_{n^{-1}}({\mathbf{T}})\cap\mathcal{N}_{<\infty}$, that has a non-empty interior, would contain an ${\mathbf{S}}_{n}\in\mathcal{N}_{\leq m-1}$. Taking $n\rightarrow\infty$ gives a sequence in $\mathcal{N}_{\leq m-1}$ with ${\mathbf{S}}_{n}\rightarrow{\mathbf{T}}$. As the set $\mathcal{N}_{\leq m-1}$ is closed (see Theorem 12) this is a contradiction to ${\mathbf{T}}\in\mathcal{N}_{m}$. Note that this proof also shows that the set $\mathcal{N}_{m}$ is an open set in the subspace topology of $\mathcal{N}_{<\infty}$ since there must exist an $n$ with $B_{n^{-1}}({\mathbf{T}})\cap\mathcal{N}_{<\infty}\subset\mathcal{N}_{m}$. ∎

The next theorem extends results that have appeared in [43, 54] for the nonnegative rank and completely characterizes all possible typical ND ranks. A typical real rank in $\otimes_{j=1}^{k}\mathbb{R}^{p_{j}}$ refers to any $r$ where the set of real, rank-$r$ tensors in $\otimes_{j=1}^{k}\mathbb{R}^{p_{j}}$ has a non-empty interior.

Theorem 9**.**

Let $r_{Re}$ be the minimum, typical, real rank for tensors in $\mathbb{R}^{p_{1}\times\cdots\times p_{k}}$. Any $r$ with $r_{Re}\leq r\leq\mathrm{maxNDrank}$ is a typical ND rank.

Proof.

We adapt the proof provided in Theorem 3.8 of [54] by replacing all instances of the the cone of nonnegative vectors $\mathbb{R}^{p_{j}}_{+}$ with the order cone $\mathcal{C}(\mathcal{P}_{j})$. This proof does not rely on any of the structure of the cone $\mathbb{R}^{p_{j}}_{+}$; the only place in this proof where the cone comes into play is in showing that the maximum rank is a typical nonnegative rank. In the ND rank case the analogous result is shown in Theorem 8. Theorem 3.8 also depends on Lemmas 3.4-3.7 in [54]. Similarly, none of these lemmas rely on the specific structure of $\mathbb{R}_{+}^{p_{i}}$ and also apply to $\mathcal{C}(\mathcal{P}_{j})$ and the ND rank. ∎

When considering rank-$r$ factorizations over the complex numbers there is only one typical rank, also referred to as the generic rank [10]. In general, there can be more than one typical, real or nonnegative rank [53]. However, for matrices, since the maximum rank of a $p_{1}\times p_{2}$ matrix is equal the minimum typical rank of $\min(p_{1},p_{2})$, this is the only typical, real or nonnegative rank. As the maximum ND rank can be larger than $\min(p_{1},p_{2})$ there can be more than one typical ND rank in the matrix setting, as is illustrated by the following corollaries.

Corollary 1**.**

If $\mathcal{C}(\mathcal{P}_{j})\subset\mathbb{R}^{p_{j}}$ has $q_{j}$ extremal rays, $j=1,2$, and $\mathcal{C}(\mathcal{P}_{2})$ is simplicial then $r$ is a typical ND rank in $\mathcal{C}(\mathcal{P}_{1})\otimes\mathcal{C}(\mathcal{P}_{2})$ if and only if $\min(p_{1},p_{2})\leq r\leq\min(q_{1},q_{2})$.

Proof.

This is a direct consequence of Theorems 6, 8, and 9; collectively these results show that any rank greater than or equal to $\min(p_{1},p_{2})$ and less than or equal to the maximum rank of $\min(q_{1},q_{2})$ is a typical rank. ∎

Corollary 2**.**

If $\mathcal{P}=\{x_{1},\ldots,x_{p}\}$ where $x_{i}\prec x_{p}$ for all $i<p$ then $r$ is a typical rank in $\mathcal{C}(\mathcal{P})\otimes\mathcal{C}(\mathcal{P})$ if $p\leq r\leq 2(p-1)$.

Proof.

The proof is the same as in Corollary 1 except that the maximum rank is now bounded below by $2(p-1)$ in Theorem 7. ∎

8 ND Rank One and Two

The relationship between rank and nonnegative rank is especially simple when the rank is either one or two: the inequality $\text{rank}({\mathbf{T}})\leq\text{rank}_{+}({\mathbf{T}})$ is satisfied with equality. In the next two results we show that this property can be extended to the nondecreasing rank.

Theorem 10**.**

If ${\mathbf{T}}\in\mathcal{C}(\times_{j=1}^{k}\mathcal{P}_{j})$ is a rank-one tensor, then ${\mathbf{T}}$ has an ND rank of one.

Proof.

By assumption ${\mathbf{T}}=\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}\neq{\mathbf{0}}$, where it can also be assumed that every ${\mathbf{v}}^{(j)}\in\mathbb{R}^{p_{j}}_{+}$ as otherwise ${\mathbf{T}}$ would have a negative entry. Suppose for a contradiction that ${\mathbf{v}}^{(i)}$ is not monotone so that there exists $x_{i}\prec y_{i}$ in $\mathcal{P}_{i}$ with $v^{(i)}_{x_{i}}>v^{(i)}_{y_{i}}$. Choosing $z_{j}\in\mathcal{P}_{j}$ with $v^{(j)}_{z_{j}}>0$ for $j\neq i$ we find that

[TABLE]

contradicting ${\mathbf{T}}\in\mathcal{C}(\times_{i=1}^{k}\mathcal{P}_{i})$. ∎

Theorem 11**.**

If ${\mathbf{T}}\in\mathcal{N}_{<\infty}$ is a rank-two matrix then ${\mathbf{T}}$ also has an ND rank of two.

Proof.

As ${\mathbf{T}}\in\mathcal{N}_{<\infty}$ we have that ${\mathbf{T}}=\sum_{i=1}^{r}{\mathbf{a}}^{(i)}\otimes{\mathbf{b}}^{(i)}$, ${\mathbf{a}}^{(i)}\in\mathcal{C}(\mathcal{P}_{1})$, ${\mathbf{b}}^{(i)}\in\mathcal{C}(\mathcal{P}_{2})$, and ${\mathbf{T}}{\mathbf{h}}=\sum_{i=1}^{r}({\mathbf{h}}^{\intercal}{\mathbf{b}}^{(i)}){\mathbf{a}}^{(i)}\in\mathcal{C}(\mathcal{P}_{1})$ for all ${\mathbf{h}}\in\mathcal{C}(\mathcal{P}_{2})^{*}$. Similarly, ${\mathbf{h}}^{\intercal}{\mathbf{T}}\in\mathcal{C}(\mathcal{P}_{2})$ for all ${\mathbf{h}}\in\mathcal{C}(\mathcal{P}_{1})^{*}$. As $\mathcal{C}(\mathcal{P}_{2})^{*}$ is full-dimensional [41] ${\mathbf{T}}(\mathcal{C}(\mathcal{P}_{2})^{*})\subseteq\mathcal{C}(\mathcal{P}_{1})$ is a two-dimensional cone. Let ${\mathbf{v}}^{(1)},{\mathbf{v}}^{(2)}$ be the extremal rays of the two-dimensional cone $\text{col}({\mathbf{T}})\cap\mathcal{C}(\mathcal{P}_{1})$. There exist ${\mathbf{w}}^{(1)},{\mathbf{w}}^{(2)}$ such that ${\mathbf{T}}={\mathbf{v}}^{(1)}\otimes{\mathbf{w}}^{(1)}+{\mathbf{v}}^{(2)}\otimes{\mathbf{w}}^{(2)}$ and it remains to show that ${\mathbf{w}}^{(i)}\in\mathcal{C}(\mathcal{P}_{2})$ for $i=1,2$. There must exist a face of $\mathcal{C}(\mathcal{P}_{1})$ that contains ${\mathbf{v}}^{(1)}$ but not ${\mathbf{v}}^{(2)}$, since otherwise $\lambda{\mathbf{v}}^{(1)}+(1-\lambda){\mathbf{v}}^{(2)}$ would be contained in $\mathcal{C}(\mathcal{P}_{1})$ for a small enough $\lambda>1$, contradicting the assumption that ${\mathbf{v}}^{(1)}$ was extremal in $\text{col}({\mathbf{T}})\cap\mathcal{C}(\mathcal{P}_{1})$. Let ${\mathbf{h}}\in\mathcal{C}(\mathcal{P}_{1})^{*}$ be a dual vector corresponding to the hyperplane that supports this face. Then ${\mathbf{h}}^{\intercal}{\mathbf{T}}=({\mathbf{h}}^{\intercal}{\mathbf{v}}^{(1)}){\mathbf{w}}^{(1)}+({\mathbf{h}}^{\intercal}{\mathbf{v}}^{(2)}){\mathbf{w}}^{(2)}=({\mathbf{h}}^{\intercal}{\mathbf{v}}^{(2)}){\mathbf{w}}^{(2)}\in\mathcal{C}(\mathcal{P}_{2})$, implying ${\mathbf{w}}^{(2)}\in\mathcal{C}(\mathcal{P}_{2})$ because ${\mathbf{h}}^{\intercal}{\mathbf{v}}^{(2)}>0$. The same argument shows that ${\mathbf{w}}^{(1)}\in\mathcal{C}(\mathcal{P}_{2})$. ∎

A ramification of the above theorems is that defining equations for the semialgebraic sets [11, Def 2.1.4] of ND rank-one tensors and ND rank-two matrices are known. In the former case, a tensor ${\mathbf{T}}$ is in $\mathcal{N}_{\leq 1}$ if and only if it is monotonic and if it satisfies the determinantal equations for a rank-one tensor [34, Sec 3.4.1]. In the latter case, a matrix ${\mathbf{T}}$ is in $\mathcal{N}_{\leq 2}$ if and only if every $3\times 3$ minor of ${\mathbf{T}}$ vanishes and if ${\mathbf{T}}$ satisfies the inequalities for the $\mathcal{H}$-representation of $\mathcal{N}_{<\infty}$ described in Sections 5 and 6. Semialgebraic conditions under which a tensor of rank two also has a nonnegative rank of two are put forward in [3]. Whenever all order cones $\mathcal{C}(\mathcal{P}_{j})$ are simplicial such conditions can be easily translated, via Theorem 4, into conditions for a tensor to have an ND rank of two.

9 ND Border Rank Equals ND Rank

Due to the presence of noise, in many applications it is unlikely that a data tensor has an exact, low-rank, ND factorization. Instead a solution to an optimization problem of the form $\text{argmin}_{\bm{\theta}\in\mathcal{N}_{\leq r}}D({\mathbf{T}},\bm{\theta})$ is sought, where $D$ is a divergence that measures the discrepancy between the observed data tensor ${\mathbf{T}}$ and an approximating, low ND rank tensor $\bm{\theta}$. To ensure that a solution to this optimization problem exists it is important that the set $\mathcal{N}_{\leq r}$ be closed. The set of tensors with real tensor rank at most $r$ is not closed, which necessitates the introduction of the concept of border rank [34, Sec 2.4.5]. It is shown in this section, extending the proof for nonnegative tensor ranks provided in [38], that $\mathcal{N}_{\leq r}$ is closed and there is no need for the additional notion of border rank.

Theorem 12** ([38] Thm 6.1).**

The set $\mathcal{N}_{\leq r}$ is closed. Equivalently, for any ${\mathbf{T}}\in\mathbb{R}^{p_{1}\times\cdots\times p_{k}}$ there exists a solution to the optimization problem $\inf_{\bm{\theta}\in\mathcal{N}_{\leq r}}\|{\mathbf{T}}-\bm{\theta}\|_{F}$, where $\|\cdot\|_{F}$ is the Frobenius norm.

Proof.

Assume that $\bm{\theta}^{(n)}=\sum_{i=1}^{r}\lambda_{i}^{(n)}\otimes_{j=1}^{k}{\mathbf{v}}_{ij}^{(n)}$ is a sequence of tensors in $\mathcal{N}_{\leq r}$ with $\|{\mathbf{T}}-\bm{\theta}^{(n)}\|_{F}\rightarrow\inf_{\bm{\theta}\in\mathcal{N}_{\leq r}}\|{\mathbf{T}}-\bm{\theta}\|_{F}$. We will first show that $\bm{\theta}^{(n)}$ has a convergent subsequence that converges to a point in $\mathcal{N}_{\leq r}$. By the continuity of the objective function this point of convergence will be a solution to the given optimization problem. The ${\mathbf{v}}_{ij}^{(n)}\in\mathcal{C}(\mathcal{P}_{j})$ are assumed, without loss of generality, to be scaled so that $\|{\mathbf{v}}_{ij}^{(n)}\|_{2}=1$ and $\lambda_{i}^{(n)}\geq 0$. It is claimed that there exists an $M$ such that $\sup_{n}\lambda_{i}^{(n)}\leq M$ for all $i$. To show this, note that $\|{\mathbf{T}}-\bm{\theta}^{(n)}\|_{F}\geq\|\bm{\theta}^{(n)}\|_{F}-\|{\mathbf{T}}\|_{F}$ and

[TABLE]

where the last inequality follows because $\langle{\mathbf{v}}_{ij}^{(n)},{\mathbf{v}}_{i^{\prime}j}^{(n)}\rangle\geq 0$ for every $i,i^{\prime},j$ as $\mathcal{C}(\mathcal{P}_{j})\subseteq\mathbb{R}_{+}^{p_{j}}$. If $\lambda_{i}^{(n)}\rightarrow\infty$ for any $i$ then $\|{\mathbf{T}}-\bm{\theta}^{(n)}\|_{F}\rightarrow\infty$, a contradiction. It follows that there exists a convergent subsequence with $\lambda_{i}^{(n_{m})}\rightarrow\lambda_{i}$ and ${\mathbf{v}}_{ij}^{(n_{m})}\rightarrow{\mathbf{v}}_{ij}\in\mathcal{C}(\mathcal{P}_{j})$, with the latter inclusion following from $\mathcal{C}(\mathcal{P}_{j})$ being closed. Thus, $\bm{\theta}^{(n_{m})}\rightarrow\tilde{\bm{\theta}}=\sum_{i=1}^{r}\lambda_{i}\otimes_{j=1}^{k}{\mathbf{v}}_{ij}\in\mathcal{N}_{\leq r}$ is a convergent subsequence. The point of convergence $\tilde{\bm{\theta}}=\text{argmin}_{\bm{\theta}\in\mathcal{N}_{\leq r}}\|{\mathbf{T}}-\bm{\theta}\|_{F}$ is a solution to the optimization problem of interest. To show that $\mathcal{N}_{\leq r}$ is closed, first observe that if ${\mathbf{T}}$ is in the closure of $\mathcal{N}_{\leq r}$, then $\inf_{\bm{\theta}\in\mathcal{N}_{\leq r}}\|{\mathbf{T}}-\bm{\theta}\|_{F}=0$. However, the point $\tilde{\theta}\in\mathcal{N}_{\leq r}$, obtained from a convergent subsequence as described above, is a solution to the minimization problem, implying that $\|{\mathbf{T}}-\tilde{\bm{\theta}}\|_{F}=\inf_{\bm{\theta}\in\mathcal{N}_{\leq r}}\|{\mathbf{T}}-\bm{\theta}\|_{F}=0$. That is, ${\mathbf{T}}=\tilde{\bm{\theta}}\in\mathcal{N}_{\leq r}$. ∎

10 Finding Low ND Rank Approximations

Solving the optimization problem $\inf_{\bm{\theta}\in\mathcal{N}_{\leq r}}\|{\mathbf{T}}-\bm{\theta}\|_{F}$ introduced in the previous section is a one possibility for finding a low ND rank approximation to a data tensor ${\mathbf{T}}$. However, when the entries of ${\mathbf{T}}$ consist of count data or positive data the Frobenius norm objective function may not be the most appropriate criteria to minimize. We examine a likelihood-based approach for estimating the mean of ${\mathbf{T}}$ in this section. It is assumed that the entries $T_{i_{1}\ldots i_{k}}$ are independently sampled from an exponential family of distributions with mean parameter [14, Ch 3] $E(T_{i_{1}\ldots i_{k}})=\theta_{i_{1}\ldots i_{k}}$ for every $(i_{1},\ldots,i_{k})\in\times_{j=1}^{k}[p_{j}]$, where $E(\cdot)$ is the expected value. The tensor $\bm{\theta}$ is constrained to have an ND rank of at most $r$ and is estimated by maximizing the observed likelihood

[TABLE]

The function $p({\mathbf{T}}|\bm{\theta})$ is the joint probability distribution or probability mass function of the tensor ${\mathbf{T}}$. Different distributional assumptions about $T_{i_{1}\ldots i_{k}}$ lead to different choices of $p(T_{i_{1}\ldots i_{k}}|\theta_{i_{1}\ldots i_{k}})$ and hence different optimization problems. Of primary interest are the following standard distributions:

If $T_{i_{1}\ldots i_{k}}\in\mathbb{R}$ then we may assume that $T_{i_{1}\ldots i_{k}}\sim\mathcal{N}(\theta_{i_{1}\ldots i_{k}},\sigma^{2})$ are independent Gaussians. The corresponding optimization problem is

[TABLE]

Note that the value of the variance $\sigma^{2}$ does not change the above optimization problem. 2. 2.

If the entries of ${\mathbf{T}}$ are counts in $\mathbb{N}$ that sum to $n$ then we may assume that ${\mathbf{T}}\sim\text{Multinomial}(\bm{\theta},n)$. Note that in this case the entries of ${\mathbf{T}}$ are not independent. The corresponding optimization problem over the probability simplex $\Delta_{\prod_{i}p_{i}-1}\subset\mathbb{R}^{p_{1}\times\cdots p_{k}}$ is

[TABLE] 3. 3.

If the $T_{i_{1}\ldots i_{k}}\in\mathbb{N}$ are counts that do not necessarily have to sum to $n$ we may assume that $T_{i_{1}\ldots i_{k}}\sim\text{Poisson}(\theta_{i_{1}\ldots i_{k}})$ independently. The corresponding optimization problem is

[TABLE]

In both this problem and the multinomial problem whenever $T_{i_{1}\ldots i_{k}}=0$ the respective $T_{i_{1}\ldots i_{k}}\log(\theta_{i_{1}\ldots i_{k}})$ term does not appear in the objective function. 4. 4.

If $T_{i_{1}\ldots i_{k}}\in(0,\infty)$ are positive then we may assume that $T_{i_{1}\ldots i_{k}}\sim\text{Exponential}(\theta_{i_{1}\ldots i_{k}}^{-1})$ independently. The corresponding optimization problem is

[TABLE]

These four optimization problems represent special cases of $\beta$-divergences, that are commonly applied to NMF problems [25], where $\beta$ is respectively equal to $2$, $1$, and [math] under the distributional assumptions in $1.$ $3.$ and $4.$ above. Similar to Proposition 7.2 in [38], which shows that optimal nonnegative approximations exist for Bregman divergences, the next result shows that this also holds for nondecreasing factorizations, but does not require any additional constraints on $\bm{\theta}$.

Lemma 10**.**

If the entries of ${\mathbf{T}}$ are in the respective subsets $\mathbb{R},\mathbb{N},\mathbb{N},(0,\infty)$ in each of the four optimization problems above then there exists a minimizer in each of the problems.

Proof.

This follows for $1.$ by Theorem 12. For the remaining three problems it is seen that there exists an $M$ where any minimizer must be contained in the region $\mathcal{S}=\{\bm{\theta}:\theta_{i_{1}\ldots i_{k}}\leq M,\forall i_{1},\ldots,i_{k}\}$, as the objective functions in $3.$ and $4.$ diverge to $\infty$ whenever $\theta_{i_{1}\ldots i_{k}}\rightarrow\infty$. In case $2.$ $\bm{\theta}$ is trivially contained in such an $\mathcal{S}$ because the probability simplex is compact. As $\mathcal{S}\cap\mathcal{N}_{\leq r}$ is compact by Theorem 12 and the objective functions are continuous on this region, existence of the solutions follows. Here the objective functions are extended continuously to take values in $\mathbb{R}\cup\{\infty\}$ with $-T_{i_{1}\ldots i_{k}}\log(0)\coloneq\infty$ and $T_{i_{1}\ldots i_{k}}/0\coloneq\infty$ when $T_{i_{1}\ldots i_{k}}\neq 0$. ∎

Before examining a least squares optimization algorithm for finding rank-$r$ ND approximations we discuss a couple of results in low-rank cases where finding an ND approximation is the same as simply performing a singular value decomposition.

The Perron-Frobenius theorem [8, Thm 3.2] that states that nonnegative, square matrices have nonnegative leading eigenvectors, has been extended to nonnegative, cubical tensors in [37, Thm 1]. The next result is reminiscent of the Perron-Frobenius theorem, but is instead concerned with the singular vectors and values of a matrix rather than the eigenvectors and eigenvalues.

Lemma 11**.**

If every fibre ${\mathbf{T}}_{i_{1}\ldots i_{j-1}\bullet i_{j+1}\ldots i_{k}}$ is in $\mathcal{C}(\mathcal{P}_{j})$ then the best rank-one approximation to a non-zero ${\mathbf{T}}$ with respect to any of the log-likelihoods mentioned in this section is equal to the best ND rank-one approximation:

[TABLE]

In particular, if ${\mathbf{T}}$ is a matrix with every column in $\mathcal{C}(\mathcal{P}_{1})$ and every row in $\mathcal{C}(\mathcal{P}_{2})$ then, up to sign changes, the first left and right singular vectors can be chosen to be in $\mathcal{C}(\mathcal{P}_{1})$ and $\mathcal{C}(\mathcal{P}_{2})$ respectively. The largest singular value of ${\mathbf{T}}$ has multiplicity one and the best ND rank-one approximation with respect to the Frobenius norm is unique.

Proof.

Throughout this proof it is assumed that the entries of ${\mathbf{T}}$ are in the sets outlined in Lemma 10 so that there exists a solution to the optimization problem.

Let $\bm{\theta}=\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}$ be a rank-one tensor, where we seek vectors ${\mathbf{v}}^{(j)}\in\mathbb{R}^{p_{j}}$ that maximize the likelihood of interest. Under the Gaussian likelihood, upon taking a gradient with respect to ${\mathbf{v}}^{(j)}$, the first-order conditions, without any monotonicity constraints, for $\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}$ to be optimal are

[TABLE]

We may assume that any optimal rank-one approximation has ${\mathbf{v}}^{(j)}\geq{\mathbf{0}}$ as the tensor $\max(\otimes_{j=1}^{k}{\mathbf{v}}^{(j)},{\mathbf{0}})=\otimes_{j=1}^{k}\max({\mathbf{v}}^{(j)},{\mathbf{0}})$, where the maximum is taken elementwise, is closer to ${\mathbf{T}}$ than $\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}$ is when $\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}$ has a negative entry. The optimal ${\mathbf{v}}^{(j)}$ is non-zero when ${\mathbf{T}}$ is non-zero since if $T_{i_{1}\ldots i_{k}}\neq 0$ then $\bm{\theta}=T_{i_{1}\ldots i_{k}}\otimes_{j=1}^{k}{\mathbf{e}}_{i_{j}}$ has a larger likelihood value than $\bm{\theta}={\mathbf{0}}$. Equation (10) shows that ${\mathbf{v}}^{(j)}$ must be a conic combination of the monotone fibres ${\mathbf{T}}_{i_{1}\ldots i_{j-1}\bullet i_{j+1}\ldots i_{k}}$, and so it must itself be in $\mathcal{C}(\mathcal{P}_{j})$.

When $k=2$ and ${\mathbf{v}}^{(1)},{\mathbf{v}}^{(2)}$ are any choice of the first left and right singular vectors, by the Eckhart-Young theorem, ${\mathbf{v}}^{(1)}\otimes{\mathbf{v}}^{(2)}$ is a best rank-one approximation of ${\mathbf{T}}$. Therefore, by the above derivation, ${\mathbf{v}}^{(j)}\in\mathcal{C}(\mathcal{P}_{j})$, up to a possible sign change. If the largest singular value was repeated then we have pairs ${\mathbf{v}}^{(j)},{\mathbf{w}}^{(j)}$ of left ($j=1$) and right ($j=2$) singular vectors. As $\mathcal{C}(\mathcal{P}_{1})$ is proper, the subspace $\text{span}({\mathbf{v}}^{(1)},{\mathbf{w}}^{(1)})$ contains a unit norm vector ${\mathbf{z}}^{(1)}\notin\pm\mathcal{C}(\mathcal{P}_{1})$. Defining the $p_{j}\times 2$ matrix ${\mathbf{R}}^{(j)}=[{\mathbf{v}}^{(j)}|{\mathbf{w}}^{(j)}]$, consider the rotation matrix ${\mathbf{V}}$ constructed so that the first column of ${\mathbf{R}}^{(1)}{\mathbf{V}}$ is equal to ${\mathbf{z}}^{(1)}$. Set ${\mathbf{z}}^{(2)}$ to be the first column of ${\mathbf{R}}^{(2)}{\mathbf{V}}$. We have that ${\mathbf{R}}^{(1)}({\mathbf{R}}^{(2)})^{\intercal}={\mathbf{R}}^{(1)}{\mathbf{V}}{\mathbf{V}}^{\intercal}({\mathbf{R}}^{(2)})^{\intercal}$ is a best rank-two approximation of ${\mathbf{T}}$ and thus ${\mathbf{z}}^{(1)}\otimes{\mathbf{z}}^{(2)}$ is a best rank-one approximation of ${\mathbf{T}}$. This contradicts the above fact that ${\mathbf{z}}^{(1)}\in\pm\mathcal{C}(\mathcal{P}_{1})$ for any best rank-one approximation.

Under the multinomial model the unconstrained maximum likelihood estimator (MLE) is well-known [1, Sec 9.6] to be equal to the averaged marginal distribution along each mode of the tensor:

[TABLE]

It is immediate that ${\mathbf{v}}^{(j)}\in\mathcal{C}(\mathcal{P}_{j})$. The Poisson MLE is equal to $c\otimes_{j=1}^{k}{\mathbf{v}}^{(j)}$ where each ${\mathbf{v}}^{(j)}$ is given as above, and $c=\langle{\mathbf{T}},{\mathbf{1}}\rangle$ is the sum of all of the entries in ${\mathbf{T}}$.

Finally, taking the gradient with respect to ${\mathbf{v}}^{(j)}$ for the exponential likelihood yields the equations

[TABLE]

This is a conic combination of ${\mathbf{T}}_{i_{1}\ldots i_{j-1}\bullet i_{j+1}\ldots i_{k}}$ and hence is monotone. Note that none of the entries of any ${\mathbf{v}}^{(l)}$ can equal zero at an optimal solution, ensuring that the above first-order condition is well-defined. This is because we define the objective function to be continuous as $\theta_{i_{1}\ldots i_{k}}$ approaches [math], and $\lim_{\theta_{i_{1}\ldots i_{k}}\rightarrow 0}\big(\log(\theta_{i_{1}\ldots i_{k}})+\tfrac{T_{i_{1}\ldots i_{k}}}{\theta_{i_{1}\ldots i_{k}}}\big)=\infty$. ∎

We remark that the constraint that each tensor fibre be in $\mathcal{C}(\mathcal{P}_{j})$ can easily be checked and is weaker than requiring that ${\mathbf{T}}\in\mathcal{N}_{<\infty}$ (Theorem 1).

The following lemma will be applied in the next section to easily find a rank-two factorization. In the case of NMF this observation was used for hierarchical clustering in [30, Sec 2.3].

Lemma 12**.**

If ${\mathbf{T}}$ has an optimal, unconstrained, rank-$r$ decomposition ${\mathbf{T}}^{(r)}=\text{argmin}_{\bm{\theta}:\mathrm{rank}(\bm{\theta})\leq r}\|{\mathbf{T}}-\bm{\theta}\|_{F}^{2}$ where ${\mathbf{T}}^{(r)}\in\mathcal{N}_{\leq r}$ then ${\mathbf{T}}^{(r)}=\mathrm{argmin}_{\bm{\theta}\in\mathcal{N}_{\leq r}}\;\|{\mathbf{T}}-\bm{\theta}\|_{F}^{2}$. For matrices with ${\mathbf{T}}^{(2)}\in\mathcal{N}_{<\infty}$ we have that ${\mathbf{T}}^{(2)}=\mathrm{argmin}_{\bm{\theta}\in\mathcal{N}_{\leq 2}}\;\|{\mathbf{T}}-\bm{\theta}\|_{F}^{2}$.

Proof.

The first statement is clear from

[TABLE]

The second statement follows from Theorem 11. ∎

There are a wide variety of algorithms for computing approximate nonnegative factorizations, the two main classes of which are multiplicative update (MU) algorithms and least squares algorithms; see [29, Ch 8] for a comprehensive account. Below we propose an algorithm for minimizing the Frobenius norm between a data tensor and a low ND rank tensor. This procedure is closely related to the hierarchical alternating least squares (HALS) algorithm [16] and is based off minimizing the following expression

[TABLE]

with respect to ${\mathbf{v}}^{(st)}\in\mathcal{C}(\mathcal{P}_{t})$, while holding all other vectors fixed. We have defined $\tilde{{\mathbf{T}}}\coloneq{\mathbf{T}}-\sum_{i\neq s}\otimes_{j=1}^{k}{\mathbf{v}}^{(ij)}$ in (11). Defining

[TABLE]

to be the unconstrained minimizer of $v_{l}^{(st)}$, the objective function in (11) can be written as

[TABLE]

where only the second term involves ${\mathbf{v}}^{(st)}$. This decomposition arises from a sum of squares decomposition in linear regression: $\|{\mathbf{y}}-{\mathbf{x}}\beta\|^{2}_{F}=\|{\mathbf{y}}-{\mathbf{x}}\hat{\beta}\|^{2}_{F}+\|{\mathbf{x}}\hat{\beta}-{\mathbf{x}}\beta\|^{2}_{F}$, with the roles of ${\mathbf{y}},{\mathbf{x}},\hat{\beta}$, and $\hat{\beta}$ respectively being played by $\tilde{{\mathbf{T}}}_{\bullet\cdots\bullet l\bullet\cdots\bullet},\otimes_{j\neq t}{\mathbf{v}}^{(sj)},\tilde{v}_{l}^{(st)}$, and $v_{l}^{(st)}$. Importantly, the second term in (12) is separable in $v_{l}^{(st)}$, $l\in[p_{t}]$, indicating that the efficient, pool-adjacent-violators algorithm (PAVA) [19, 57] can be used to minimize this quadratic program over $\mathcal{C}(\mathcal{P}_{t})$. Our ND hierarchical least squares algorithm is a block-coordinate descent routine that proceeds by cycling through the PAVA updates for each vector ${\mathbf{v}}^{(ij)}$ until convergence.

Due to the Frobenius norm optimization problem being jointly non-convex in the ${\mathbf{v}}^{(ij)}$s there is no guarantee that the sequence of updates converges to a global optimum. However, Proposition 2.7.1 of [9] ensures that any limit point of the ${\mathbf{v}}^{(ij)}$ iterates is a stationary point as long as none of the vectors ${\mathbf{v}}^{(ij)}$ are updated to the zero vector over the course of the algorithm. The issue that can occur when any ${\mathbf{v}}^{(sj)}={\mathbf{0}}$ is that the componentwise objective function (12) is not strictly convex because the relevant quadratic term vanishes. As there can exist multiple stationary points, it is recommended that multiple, distinct initializations be chosen in the ND HALS algorithm. One possible choice of initialization is to first run an unconstrained, rank-$r$ approximation algorithm on ${\mathbf{T}}$, such as alternating least squares [36], to obtain vectors $\tilde{{\mathbf{v}}}^{(ij)}$. These vectors can then respectively be orthogonally projected onto the cone $\mathcal{C}(\mathcal{P}_{j})$ by again using PAVA to solve ${\mathbf{v}}^{(ij)}=\text{argmin}_{{\mathbf{v}}^{(ij)}\in\mathcal{C}(\mathcal{P}_{j})}\|\tilde{{\mathbf{v}}}^{(ij)}-{\mathbf{v}}^{(ij)}\|_{2}^{2}$.

11 Applications of Low ND Rank Decompositions

In this section two applications of low ND rank factorizations are presented.

11.1 Pig Weight Data

The pig weight dataset [20, Sec 3.2] available in the fds R package [48] records the weights of $48$ growing pigs over the course of nine weeks. Figure 5(a) illustrates that all of the pigs have monotonically increasing weights over time. If ${\mathbf{T}}\in\mathbb{R}^{48\times 9}$ is the matrix of weights we may assume that the poset corresponding to the columns is a chain, while no ordering constraints, apart from nonnegativity, are placed on the rows. If ${\mathbf{T}}=\sum_{i=1}^{r}{\mathbf{a}}^{(i)}({\mathbf{b}}^{(i)})^{\intercal}$ is a rank-$r$ ND factorization, the weight profile of pig $j$ is ${\mathbf{T}}_{j\cdot}=\sum_{i=1}^{r}a^{(i)}_{j}({\mathbf{b}}^{(i)})^{\intercal}$. Up to scale factors, the vectors ${\mathbf{b}}^{(i)}$ can be interpreted as growth curves for $r$ different “subpopulations” of pigs, while $a^{(i)}_{j}$ indicates the degree of membership of pig $j$ to subpopulation $i$. The growth curve of each pig is modeled to approximately be a conical combination of the subpopulation curves. If only an NMF was applied to this dataset the ${\mathbf{b}}^{(i)}$ may not be as easily interpretable as growth curves. Furthermore, the second, right, singular vector certainly cannot be interpreted as a growth curve as it does not even contain entries that have the same sign.

To fit this dataset we find an ND rank two approximation ${\mathbf{T}}\approx{\mathbf{a}}^{(1)}\otimes{\mathbf{b}}^{(1)}+{\mathbf{a}}^{(2)}\otimes{\mathbf{b}}^{(2)}$ with ${\mathbf{a}}^{(i)}\in\mathbb{R}^{48}_{+}$ and ${\mathbf{b}}^{(i)}\in\mathcal{C}(\{1,\ldots,9\})$. A rank of two was chosen because this rank is sufficient to yield a low-reconstruction error for ${\mathbf{T}}$; the first four squared singular values of ${\mathbf{T}}$ are $121{,}000$, $613$, $129$, and $109$. The rank-two SVD approximation ${\mathbf{T}}^{(2)}$ of ${\mathbf{T}}$ has nonnegative entries and monotone increasing rows, implying that ${\mathbf{T}}^{(2)}\in\mathcal{N}_{\leq 2}$ (Theorem 11). The vectors ${\mathbf{b}}_{1},{\mathbf{b}}_{2}$ can be taken to be the pair of rows of ${\mathbf{T}}^{(2)}$ that have the largest angle between each other. Such a pair can easily be found by computing dot products between the rows of ${\mathbf{T}}^{(2)}$ that are rescaled to have unit norm. This approach has been used in the NMF literature for matrices that have a nonnegative rank of two which have the property of being separable [29, Sec 4.1]. The corresponding coefficients ${\mathbf{a}}^{(1)},{\mathbf{a}}^{(2)}$ are found by a linear regression. Lastly, we rescale ${\mathbf{b}}^{(1)},{\mathbf{b}}^{(2)}$ so that the mean weight of each ${\mathbf{b}}^{(i)}$ is equal to the average weight across all the pigs: $\tfrac{1}{9}\sum_{l=1}^{9}b^{(i)}_{l}=\tfrac{1}{9*48}\sum_{i=1}^{48}\sum_{j=1}^{9}T_{ij}$. Figure 5(a) displays the two estimated growth curves ${\mathbf{b}}^{(1)}$ and ${\mathbf{b}}^{(2)}$, overlaid on the observed growth curves. Interestingly, around week five the slope of both estimated curves change, a phenomena that may warrant further study. In Figure 5(b) we see that the fitted curve $a_{i}^{(1)}{\mathbf{b}}^{(1)}+a_{i}^{(2)}{\mathbf{b}}^{(2)}$, for pig $i$, is nearly a convex combination of ${\mathbf{b}}^{(1)}$ and ${\mathbf{b}}^{(2)}$, as the scale factor $a^{(1)}_{i}+a^{(2)}_{i}$ is close to one. The relative importance of ${\mathbf{b}}^{(1)}$ as compared to ${\mathbf{b}}^{(2)}$ noticeably varies across the pigs.

One issue that has not yet been discussed in this work is that ND factorizations are not unique in general, even after accounting for scaling and permutation of the rank-one factors. Various regularization penalties or additional constraints have been introduced in the NMF literature [29, Ch 4] to select a factorization that has desirable properties. In the above analysis we have implicitly selected the minimum volume [40] ND factorization by choosing rows from ${\mathbf{T}}^{(2)}$ to serve as the ${\mathbf{b}}^{(i)}$. An alternative here is take a maximum volume factorization by choosing the ${\mathbf{b}}^{(i)}$ to be the extremal rays of the cone $\mathcal{C}(\{1,\ldots,9\})\cap\text{row}({\mathbf{T}}^{(2)})$. Either factorization will yield the same reconstruction error $\|{\mathbf{T}}-{\mathbf{T}}^{(2)}\|^{2}_{F}$, but will lead to potentially different interpretations with respect to the ${\mathbf{b}}^{(i)}$ curves and the ${\mathbf{a}}^{(i)}$ membership weights.

11.2 Perceived Mental Health of Canadians

We examine a subset of the data from the annual Canadian Community Health Survey (CCHS) [52] from $2016$-$2022$ that records responses of whether an individual perceives their mental health to be either fair or poor. The survey also records the sex and age group of the participant. The $(5,7,2)$-dimensional tensor ${\mathbf{T}}$ in Table 2 displays the proportion of respondents who indicated that they had either poor or fair mental health.

Using the ND HALs algorithm from Section 10 we fit both rank-two and rank-one approximations to the data tensor of the form $\sum_{i=1}^{r}{\mathbf{a}}^{(i)}\otimes{\mathbf{b}}^{(i)}\otimes{\mathbf{c}}^{(i)}$. The residual sum of squares for the respective approximations are $145$ and $55$ with the total sum of squares equaling $6925$. We only impose nonnegativity constraints on the sex mode, while the Hasse diagrams in Figure 6 provides orderings on the age group and years that are consistent with the observed data. Specifically, the mental health of respondents deteriorated in 2020 and the subsequent two years, conceivably as a result of the COVID-19 pandemic. Respondents in the second youngest age group were most likely to perceive their mental health as fair or poor while respondents in the oldest age group were most likely to not provide this response.

There is a marked difference between how individuals in the 12-17 age group responded between females and males, with females having poorer perceived mental health. In fitting a rank-two ND decomposition the two rank-one factors provide a low-dimensional summary of this sex-age interaction. The rank-one factors of the two rank-one terms are shown in Table 3. These factors are scaled so that $\|{\mathbf{c}}^{(i)}\|_{1}=1$ and $\|{\mathbf{b}}^{(i)}\|_{1}=7$. We see that the factor ${\mathbf{a}}^{(2)}\otimes{\mathbf{b}}^{(2)}\otimes{\mathbf{c}}^{(2)}$ only alters the female portion of the tensor and represents a deviation between the male and female responses. In comparison with ${\mathbf{a}}^{(1)}$, the 12-17 age group in ${\mathbf{a}}^{(2)}$ takes on larger values relative to the 18-34 age group. Furthermore, the impact of pandemic is more apparent in ${\mathbf{b}}^{(2)}$ as compared to ${\mathbf{b}}^{(1)}$. In summary, the low ND rank decomposition allows for a concise interpretation of the data tensor and illustrates that young females in 2022 were particularly prone to respond that they had fair or poor mental health as compared to young males.

12 Conclusion and Future Directions

In this article the concept of the nondecreasing rank was introduced and properties of nondecreasing factorizations were examined. There remains a wealth of relevant, open questions to be explored, a few of which we mention here. One property of ND factorizations not extensively discussed is the uniqueness of such factorizations. Like with nonnegative matrix factorizations, it is expected that only under similarly strong conditions, such as separability [29, Sec 4.2.2], will a matrix have a unique ND factorization. Also not discussed, is how to choose a suitable rank $r$ for performing an ND factorization. An easy method for doing this is to plot the value of a divergence objective function against $r$ and look for a kink in the plot where the divergence stops decreasing rapidly in response to an increase of $r$ [15]. In Section 10 an algorithm is introduced for minimizing the Frobenius norm between a tensor and its ND approximation. Developing efficient multiplicative update [35] algorithms for finding ND approximations would be a valuable contribution. To construct such an algorithm more constraints as compared to NMF, stemming from the monotonicity constraints of the order cones, have to be contended with.

Nondecreasing factorizations are closely related to ideas in the field of order constrained statistical inference [7, 47]. Beyond monotonicity constraints, related shape constraints for functions, such as convexity and unimodality, would be interesting to inspect for low-rank structures. Unimodal matrix factorizations have previously been examined in [4, 13]. Unimodality is related to the umbrella order where $x_{1}\prec\cdots\prec x_{l-1}\prec x_{l}\succ x_{l+1}\cdots\succ x_{p}$. The main distinction is that in a unimodal vector it is not known which of the entries is the maximum, while for the umbrella order above the $x_{l}$th entry is posited to be the maximum. The non-convex cone of unimodal vectors is therefore equal to a union of umbrella-order cones. Problems involving convex, unimodal, and monotone factorizations can also be contextualized from an infinite-dimensional, functional data analysis perspective, where notions of convex sets and rank continue to be well-defined.

Acknowledgments

The author acknowledges the support of the Natural Sciences and Engineering Research Council of Canada (NSERC), [RGPIN-2025-03968, DGECR-2025-00237].

Cette recherche a été financée par le Conseil de recherches en sciences naturelles et en génie du Canada (CRSNG), [RGPIN-2025-03968, DGECR-2025-00237].