Tom Szwagier
email: tom.szwagier@inria.fr

Xavier Pennec
email: xavier.pennec@inria.fr

Rethinking the Riemannian Logarithm on Flag Manifolds as an Orthogonal Alignment Problem

Keywords: Flag manifolds, Riemannian logarithm, Orthogonal alignment, Procrustes analysis

Flags are sequences of nested linear subspaces of increasing dimension. They belong to smooth manifolds generalizing Grassmannians and bring a richer multi-scale point of view to the traditional subspace methods in statistical analysis. Hence, there is an increasing interest in generalizing the formulae and statistical methods already developed for Grassmannians to flag manifolds. In particular, it is critical to compute accurately and efficiently the geodesic distance and the logarithm due to their fundamental importance in geometric statistics. However, there is no explicit expression known in the case of flags. In this work, we exploit the homogeneous quotient space structure of flag manifolds and rethink the geodesic endpoint problem as an alignment of orthogonal matrices on their equivalence classes. The relaxed problem with the Frobenius metric surprisingly enjoys an explicit solution. This is the key to modify a previously proposed algorithm. We show that our explicit alignment step brings drastic improvements in accuracy, speed and radius of convergence, in addition to overcoming the combinatorial issues raised by the non-connectedness of the equivalence classes.

Introduction

Flags are sequences of nested linear subspaces of increasing dimension. They are important in statistical analysis [START_REF] Ma | The Flag Manifold as a Tool for Analyzing and Comparing Sets of Data Sets[END_REF][START_REF] Pennec | Barycentric Subspace Analysis on Manifolds[END_REF] due to the multi-scale information they provide, compared to traditional subspace methods involving Grassmann manifolds [START_REF] Bendokat | A Grassmann Manifold Handbook: Basic Geometry and Computational Aspects[END_REF]. Flags of a given type form a Riemannian manifold that generalizes Grassmannians. Hence there is a natural interest in generalizing the formulae and statistical methods already developed on Grassmannians.

The geodesic distance and logarithm are central tools in statistics on Riemannian manifolds, as they allow notably to discriminate, interpolate, and optimize [START_REF] Pennec | Riemannian Geometric Statistics in Medical Image Analysis[END_REF]. Their explicit formulae are known for Grassmannians, and are related to the problem of finding principal vectors and angles between linear subspaces, which can be solved using Singular Value Decomposition (SVD) [START_REF] Bendokat | A Grassmann Manifold Handbook: Basic Geometry and Computational Aspects[END_REF]. However no explicit formula is known for flags. One of the main hurdles seems to be the relatively complex structure of the tangent space, that must carry the information of nestedness of several subspaces, compared to Grassmannians which account for only two subspaces, and yield closed-form expressions for the geodesics in terms of sines and cosines of principal angles, like in [START_REF] Bendokat | A Grassmann Manifold Handbook: Basic Geometry and Computational Aspects[END_REF]Prop. 3.3 ].

The question of the distance on flag manifolds is first addressed in [12, Prop. 10 ], but the proposition is actually misleading. Indeed, the proof of this proposition is based on the implicit knowledge of the horizontal tangent vector allowing to shoot from one flag to the other with minimal distance, that is the Riemannian logarithm. However, if this vector is known in the geodesic shooting problem, it is actually unknown in the geodesic endpoint problem, otherwise the geodesic distance would also be known, as it is the norm of the logarithm. Without clearly stating this assumption, the proposition is false. For instance, two orthogonal matrices representing the same flag would have a positive distance. Subsequent works rather focus on algorithms to approximate the Riemannian logarithm and distance, either based on a quotient space point of view and using alternated projections on the horizontal and vertical tangent spaces [START_REF] Ma | The Flag Manifold as a Tool for Analyzing and Comparing Sets of Data Sets[END_REF][START_REF] Ma | Self-organizing mappings on the flag manifold with applications to hyper-spectral image data analysis[END_REF], or following an optimization approach [START_REF] Nguyen | Closed-form Geodesics and Optimization for Riemannian Logarithms of Stiefel and Flag Manifolds[END_REF][START_REF] Nijhawan | Flag Manifold-Based Precoder Interpolation Techniques for MIMO-OFDM Systems[END_REF].

In this work, we exploit the homogeneous quotient space structure of flag manifolds and rethink the geodesic endpoint problem (i.e., the search for a geodesic of minimal length joining two flags) as an alignment problem on the orthogonal group, that is the search for orthogonal matrices that are the closest among their equivalence classes. As it does not have an explicit solution, we relax it using the Frobenius metric. Our key result is that this constrained orthogonal Procrustes problem actually enjoys an explicit solution. We show that modifying a previously existing algorithm [4, Alg. 1 ] with this new result drastically improves its accuracy, speed and radius of convergence, in addition to overcoming the combinatorial issues raised by the non-connectedness of the stabilizer subgroup.

Flag Manifolds

In this section, we briefly define flag manifolds and the important properties for the next sections. A more complete introduction is given in [START_REF] Ye | Optimization on flag manifolds[END_REF].

Definition 1 (Flag, Signature). Let n ≥ 2 and d 1 < d 2 < • • • < d r = n be a sequence of strictly increasing natural numbers. A flag of signature (d 1 , . . . , d r ) is a sequence of properly nested linear subspaces of R n , E 1 ⊂ • • • ⊂ E r = R n of respective dimension d 1 < • • • < d r = n. Flags with signatures of length n, i.e. (1, 2, . . . , n -1, n) are called complete.
Flags can equivalently be defined as a sequence of incremental orthogonal subspaces

V 1 ⊥ • • • ⊥ V r ,
by taking the orthogonal complement of one nested subspace into the next one. One then has

E i = i j=1 V j (i = 1 . . r). The asso- ciated sequence of increments I := (n 1 , . . . , n r ) := (d 1 , d 2 -d 1 , . . . , d r -d r-1 )
is called the type of the flag. We will use this definition of flags in the rest of the paper.

⊂ ⊂

(a) as a sequence of nested subspaces A flag can be represented as an orthogonal matrix whose columns form an orthonormal basis that is adapted to the sequence of incremental orthogonal subspaces. More specifically, for a given sequence of subspaces

• Q 1 Q 2 Q 3 (b) as an orthogonal matrix
V 1 ⊥ • • • ⊥ V r such that r i=1 V i = R n , one can find an orthonormal basis Q i := q 1 i | . . . |q ni i ∈ R n×ni for each V i ,

and the concatenation of those bases forms an orthogonal matrix

Q = [Q 1 | . . . |Q r ] ∈ O(n) that is a representative of the flag V 1 ⊥ • • • ⊥ V r . A flag of type (1, 1, 1
) is represented in Figure 1, both in terms of nested subspaces and in the orthogonal representation.

Theorem 1 (Flag Manifolds [START_REF] Monk | The Geometry of Flag Manifolds[END_REF]). The set of all flags of type I is a connected compact smooth manifold. It is noted Flag (n 1 , . . . , n r ), or Flag(I) for short.

There exist several orthonormal bases adapted to a flag, hence flags are actually equivalence classes of orthogonal matrices. More precisely, let Q ∈ O(n) be a representative of a given flag and

R = diag (R 1 , . . . , R r ) ∈ O(n), with R i ∈ O(n i ) (i = 1 . . r). Then Q and QR = [Q 1 R 1 | . . . |Q r R r ]
represent the same flag. Indeed, the right multiplication by the block diagonal orthogonal matrix R on Q only rotates and reflects the incremental orthonormal bases within their subspaces, but do not change their span. Therefore, flag manifolds are homogeneous quotient spaces of O(n) [12, Eq. ( 8)]

Flag(I) ∼ = O(n)/ (O(n 1 ) × • • • × O(n r )) . (1)
For the sake of readability, we will thereafter write O(I)

:= O(n 1 ) × • • • × O(n r ). Therefore one has Flag(I) ∼ = O(n)/O(I) and R ∈ O(I).
Hence, a flag is an equivalence class of orthogonal matrices

[Q] := {QR, R ∈ O(I)} := Q • O(I) (2)
and one defines the canonical projection

Π c : Q ∈ O(n) → [Q] ∈ O(n)/O(I).
One can also embed flags in a product of Grassmannians Gr(n 1 , n) × • • • × Gr(n r , n) := Gr(I), representing them as a sequence of orthogonal projection matrices onto the sequence of incremental orthogonal subspaces of the flag. The embedding map Π Gr(I)

: [Q 1 | . . . |Q r ] ∈ O(n) → Q 1 Q 1 , . . . , Q r Q r ∈ Gr(I)
then removes the necessity to work with equivalence classes.

The theory of Riemannian submersions and quotient spaces [START_REF] O'neill | The fundamental equations of a submersion[END_REF] tells that the Riemannian manifold structure of flags can be deduced from the one of O(n) and O(I) with their canonical metric g Q , given by a multiple of the Frobenius metric: [START_REF] Guigui | Introduction to Riemannian Geometry and Geometric Statistics: From Basic Theory to Implementation with Geomstats[END_REF][START_REF] Ye | Optimization on flag manifolds[END_REF]. Their respective tangent spaces at the identity are o(n) := Skew n , the set of n × n skew-symmetric matrices, and o(I) := Skew I := diag (Skew n1 , . . . , Skew nr ). In the fiber bundle vocabulary, o(I) is referred to as the vertical space (noted Ver(I)), and its orthogonal complement in Skew n as the horizontal space, containing skew-symmetric matrices with diagonal zero blocks

g Q (A, B) = 1 2 A, B F := 1 2 tr A B [
Hor(I) :=          H :=      0 n1 H 1,2 . . . H 1,r -H 1,2 0 n2 . . . H 2,r . . . . . . . . . . . . -H 1,r -H 2,r . . . 0 nr      , H ij ∈ R ni×nj          . ( 3
)
The tangent space of flag manifolds is then given by

T [Q] (O(n)/O(I)) = Q Hor(I). (4)
The geodesics in the base space O(n)/O(I) with the canonical metric are inherited from the horizontal geodesics in the top space O(n) [START_REF] Ye | Optimization on flag manifolds[END_REF] exp

[Q] (A) = Q exp( Q A ) .
A common approach for computing the Riemannian logarithm and distance is through solving the geodesic endpoint problem arg min

H∈Hor(I),R∈O(I)

1 2 H 2 F subject to P exp(H) = QR. (5)
However, it has a priori no solution for flag manifolds. In the literature, the geodesic endpoint problems are generally optimized using a gradient descent on the tangent space [START_REF] Nguyen | Closed-form Geodesics and Optimization for Riemannian Logarithms of Stiefel and Flag Manifolds[END_REF][START_REF] Zimmermann | Computing the Riemannian Logarithm on the Stiefel Manifold: Metrics, Methods, and Performance[END_REF]. In the following section, we propose a different approach based on the previously described quotient space structure, and reformulate the geodesic endpoint problem as an alignment on the equivalence classes.

Alignment on Flag Manifolds

The notion of alignment in Riemannian quotient spaces is mentioned in Moreover, [START_REF] Ye | Optimization on flag manifolds[END_REF]Thm. 4 ] ensures that such an alignment always exists on flag manifolds. Hence we get an approach to the computation of the Riemannian logarithm that leverages more information than classical optimization algorithms solving the geodesic endpoint problem, because it takes into consideration the homogeneous quotient space structure of flag manifolds, and takes advantage of the results that already exist on the top space O(n). Notably, even if the proposed algorithm does not converge, one still has a way to interpolate between two flags using the logarithm in O(n). Just like the geodesic endpoint problem, the alignment problem [START_REF] Nguyen | Closed-form Geodesics and Optimization for Riemannian Logarithms of Stiefel and Flag Manifolds[END_REF] does not have either an analytic solution, to the best of our knowledge. However, the embedding of O(n) in R n×n raises the natural idea of relaxing the problem and working with the Frobenius metric.

P -QR 2 F . (7)
This name refers to the celebrated orthogonal Procrustes problem [START_REF] Schönemann | A generalized solution of the orthogonal procrustes problem[END_REF], which also involves the minimization of a Frobenius distance on an orthogonal group.

Here, our optimization problem is constrained to the isotropy subgroup O(I) of flag manifolds in O(n). We now give the main result of the paper.

Theorem 2 (O(I) Procrustes Solution on Flag Manifolds).

Let P := [P 1 | . . . |P r ], Q := [Q 1 | . . . |Q r ] ∈ O(n)
and let us write the SVDs

P i Q i := U i Σ i V i (i = 1 . . r), with U i , V i ∈ O(n i ) and Σ i ∈ diag R ni ≥0 . A solution R * = diag (R * 1 , . . . , R * r ) ∈ O(I) to the O(I)-constrained Procrustes problem on flag manifolds (7) is given by R * i = U i V i (i = 1 . . r) . ( 8
)
The uniqueness is conditioned on the uniqueness of the SVD.

Proof. We can show that this O(I) Procrustes problem is equivalent to the independent resolution of classical orthogonal Procrustes problem [START_REF] Schönemann | A generalized solution of the orthogonal procrustes problem[END_REF], i.e. min R∈O(I) evaluation process is as follows. First we generate an orthogonal matrix P drawn from a uniform distribution on O(n). Second we generate H ∈ Hor(I) drawn from a uniform distribution (with a specified norm 1 2 H F ∈ π 5 , π 2 , π ) and get the aligned endpoint Q aligned = P exp(H). Third we generate R ∈ O(I) drawn from a uniform distribution and get the endpoint Q = P exp(H)R. We run both algorithms with a maximal number of iterations of 10, and repeat the experiment 10 times independently, each time outputting an optimal horizontal vector H * . The evaluation metrics are: endpoint error ( Π Gr(I) (P exp(H

P -QR 2 F = r i=1 min Ri∈O(ni) P i -Q i R i 2 F .
* ))-Π Gr(I) (Q) F ), distance error ( 1 2 | H * F -H F |), horizontal error ( 1 2 H * -H F )
, number of iterations and computing time. The results are reported in Table 1 and the associated learning curves are illustrated in Figure 3.

It is clear that the introduction of O(I) Procrustes alignment as an initialization and a substitute for the vertical projection drastically improves [4, Alg. 1 ], in terms of accuracy, speed and radius of convergence, on all the evaluation metrics, as well as overcomes the difficulty of non-connectedness of the orbits. The difference of computing time is particularly important, because Algorithm 1 does not require to try all the connected components of O(I) and converges faster.

Conclusion and Perspectives

In this work, we first gave a quick overview on the state of research in flag manifolds logarithm and distance. Then, rethinking the geodesic endpoint problem as an alignment problem, we proposed a relaxed solution involving orthogonal Procrustes analysis. Finally, we showed how introducing it into [4, Alg. 1 ] drastically improves its accuracy, speed and radius of convergence, as well as overcomes the combinatorial drawback of working on fully oriented flag manifolds. The code is available on GitHub1 .

In future work, we could first improve the computational efficiency of the algorithm, by using linear algebra tricks like the 2k Embedding of [3, Thm. 1 ]. Second we could investigate the convergence properties of our algorithm, by understanding the link between the alignment of the target point on the target orbit and the gradient of the geodesic endpoint problem.

Fig. 1 :

 1 Fig. 1: Different representations of a flag of type (1, 1, 1).

[ 2 ,Definition 2 ( 2 F 6 )

 2226 Def. 5.1.6 ]. It refers to the problem of finding a pair of points which minimizes the top space distance within their respective equivalence classes. Alignment on Flag Manifolds). Let P, Q ∈ O(n) and d O(n) (P, Q) 2 := 1 2 log(P Q) , the geodesic distance on O(n). We define the alignment problem on flag manifolds as arg min R∈O(I) d O(n) (P, QR) 2 . (Given two flags [P ], [Q] ∈ O(n)/O(I), finding a point Q * := QR * which minimizes the O(n) geodesic distance to P in the orbit [Q] yields both the geodesic distance between the two flags d([P ], [Q]) = d O(n) (P, Q * ) and the logarithm log [P ] ([Q]) = P log(P Q * ), according to [2, Prop. 5.1.3 & Def. 5.1.6 ].

3. 1

 1 Relaxing Orthogonal to Frobenius Distance Definition 3 (O(I) Procrustes Problem on Flag Manifolds). Let P, Q ∈ O(n). We define the O(I)-constrained Procrustes problem on flag manifolds as arg min R∈O(I)

Fig. 3 :

 3 Fig.3: Evolution of the endpoint errors along the iterations for different shooting vector lengths1 2 H F ∈ π 5 , π 2 , π . Solid: Algorithm 1, Dashed:[START_REF] Ma | Self-organizing mappings on the flag manifold with applications to hyper-spectral image data analysis[END_REF].

Table 1 :

 1 Evaluation of Algorithm 1 and [4] on Flag(1, 2, 2, 15, 80) averaged over 10 experiments, for different geodesic distances (from top to bottom: π

	5 , π 2 , π).

Link for the code: https://github.com/tomszwagier/flag-manifold-distance.

Acknowledgements: This work was supported by the ERC grant #786854 G-Statistics from the European Research Council under the European Union's Horizon 2020 research and innovation program and by the French government through the 3IA Côte d'Azur Investments ANR-19-P3IA-0002 managed by the National Research Agency.

Algorithm 1 Flag Alignment

Fig. 2: The flag alignment Algorithm 1 and a conceptual visualization. The point Q moves in its equivalence class [Q] at each iteration, with the aim that the O(n) logarithm P X t becomes horizontal.

Introducing O(I) Procrustes into an Existing Algorithm

An algorithm to compute the Riemannian logarithm and distance on flag manifolds is proposed in [4, Alg. 1 ]. This algorithm takes as inputs two special orthogonal matrices P, Q ∈ SO(n). It aims at finding a horizontal vector H ∈ Hor(I) and a vertical vector V ∈ Ver(I) such that one can write Q = P exp(H) exp(V ).

It is based on iterated alternating projections on the horizontal and vertical spaces. One important drawback of this method is that it implicitly assumes that the endpoint Q is in the same connected component in O(I) as the horizontally aligned point Q aligned = P exp(H), which is generally not the case. The authors highlight this drawback and explain that the algorithm is actually working on fully oriented flag manifolds SO(n)/ (SO(n 1 ) × • • • × SO(n r )). To overcome this issue, they first create 2 r-1 "equivalents" of the endpoint Q in all the connected components, then run 2 r-1 times their algorithm, and finally take the best outcome among the different runs.

In this work, we take advantage of our new result (Theorem 2), which gives a global minimum on the whole isotropy group O(I), to overcome this drawback. We reinterpret [4, Alg. 1 ] as an alignment problem on flag manifolds and introduce O(I) Procrustes alignment both in the initialization and as a substitute for the vertical projection step onto o(I). The algorithm and a conceptual visualization of it are given in Figure 2.

Numerical Experiments

We now evaluate and compare Algorithm 1 to the alternating projections algorithm proposed in [START_REF] Ma | Self-organizing mappings on the flag manifold with applications to hyper-spectral image data analysis[END_REF], run on all the connected components of O(I). The general


Lek-Heng Lim

Pierre Comon
email: pcomon@unice.fr

Richard Allan

NONNEGATIVE APPROXIMATIONS OF NONNEGATIVE TENSORS

come

Introduction

The decomposition of a tensor into a minimal sum of outer products of vectors was first studied by Hitchcock [START_REF] Hitchcock | The expression of a tensor or a polyadic as a sum of products[END_REF][START_REF] Hitchcock | Multiple invariants and generalized rank of a p-way matrix or tensor[END_REF] in 1927. The topic has a long and illustrious history in algebraic computational complexity theory (cf. [START_REF] Bürgisser | Algebraic complexity theory[END_REF] and the nearly 600 references in its bibliography) dating back to Strassen's celebrated result [START_REF] Strassen | Gaussian elimination is not optimal[END_REF]. It has also recently found renewed interests, coming most notably from algebraic statistics and quantum computing.

However the study of the corresponding approximation problem, i.e. the approximation of a tensor by a sum of outer products of vectors, probably first surfaced as data analytic models in psychometrics in the work of Harshman [START_REF] Harshman | Foundations of the parafac procedure: models and conditions for an explanatory multi-modal factor analysis[END_REF], who called his model parafac (for Parallel Factor Analysis), and the work of Carrol and Chang [START_REF] Carroll | Analysis of individual differences in multidimensional scaling via n-way generalization of Eckart-Young decomposition[END_REF], who called their model candecomp (for Canonical Decomposition).

The candecomp/parafac model, sometimes abbreviated as cp model, essentially asks for a solution to the following problem: given a tensor A ∈ R d 1 ו••×d k , find an optimal rank-r approximation to A, [START_REF] Bini | O(n 2.7799 ) complexity for n × n approximate matrix multiplication[END_REF] X r ∈ argmin rank(X)≤r A -X , Key words and phrases. Nonnegative tensors, nonnegative hypermatrices, nonnegative tensor decompositions, nonnegative tensor rank, probabilistic latent semantic indexing, candecomp, parafac, tensor norm, tensor Brègman divergence.

or, more precisely, find scalars λ p and unit vectors 1 u p , v p , . . . , z p , p = 1, . . . , r, that minimizes [START_REF] Bini | Approximate solutions for the bilinear form computational problem[END_REF] A -

r p=1 λ p u p ⊗ v p ⊗ • • • ⊗ z p .
The norm • here is arbitrary and we will discuss several natural choices in the next section. When k = 2, A becomes a matrix and a solution to the problem when • is unitarily invariant is given by the celebrated Eckart-Young theorem: X r may be taken to be

X r = r p=1 σ p u p ⊗ v p ,
where σ 1 ≥ • • • ≥ σ r are the first r singular values of A and u p , v p the corresponding left and right singular vectors. However when k ≥ 3 the problem becomes more subtle. In fact, a global minimizer of (2) may not even exist as soon as k ≥ 3; in which case the problem in ( 1) is ill-posed because the set of minimizers is empty. We refer the reader to Section 5 for examples and discussions. Nevertheless we will show that for nonnegative tensors the problem of finding a best nonnegative rank-r approximation always has a solution, i.e. (2) will always have a global minimum when A and u p , v p , . . . , z p are required to be nonnegative. Such nonnegativity arises naturally in applications. For example, in the context of chemometrics, sample concentration and spectral intensity often cannot assume negative values [START_REF] Bro | A fast non-negativity constrained least squares algorithm[END_REF][START_REF] Bro | Least squares algorithms under unimodality and non-negativity constraints[END_REF][START_REF] Carroll | Fitting of the latent class model via iteratively reweighted least squares candecomp with nonnegativity constraints[END_REF][START_REF] Krijnen | Contrastvrije oplossingen van het candecomp/parafacmodel[END_REF][START_REF] Paatero | A weighted non-negative least squares algorithm for three-way parafac factor analysis[END_REF][START_REF] Paatero | Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values[END_REF]. Nonnegativity can also be motivated by the data analytic tenet [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF] that the way 'basis functions' combine to build 'target objects' is an exclusively additive process and should not involve any cancellations between the basis functions. For k = 2, this is the motivation behind nonnegative matrix factorization (nmf) [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF][START_REF] Paatero | Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values[END_REF], essentially a decomposition of a nonnegative matrix A ∈ R m×n into a sum of outer-products of nonnegative vectors,

A = W H ⊤ = r p=1 w p ⊗ h p ,
or, in the noisy situation, the approximation of a nonnegative matrix by such a sum:

min W ≥0,H≥0 A -W H ⊤ = min wp≥0,hp≥0 A - r p=1 w p ⊗ h p .
The generalization of nmf to tensors of higher order yields a model known as nonnegative parafac [START_REF] Carroll | Fitting of the latent class model via iteratively reweighted least squares candecomp with nonnegativity constraints[END_REF][START_REF] Krijnen | Contrastvrije oplossingen van het candecomp/parafacmodel[END_REF][START_REF] Paatero | A weighted non-negative least squares algorithm for three-way parafac factor analysis[END_REF], which has also been studied more recently under the name nonnegative tensor factorization (ntf) [START_REF] Shashua | Non-negative tensor factorization with applications to statistics and computer vision[END_REF]. As we have just mentioned, a general tensor can fail to have a best low-rank approximation. So the first question that one should ask in a multilinear generalization of a bilinear model is whether the generalized problem would still have a solution -and this was the question that Harshman posed. More generally, we will show that nonnegative parafac always has a solution for any continuous measure of proximity satisfying some mild conditions, e.g. norms or Brègman divergences. These include the sum-of-squares loss and Kullback-Leibler divergence commonly used in nmf and ntf.

The following will be proved in Sections 6 and 7. Let

Ω 0 ⊆ Ω ⊆ R d 1 ו••×d k + be closed convex subsets. Let d : Ω × Ω 0 →
R be a norm or a Brègman divergence. For any nonnegative tensor A ∈ Ω and any given r ∈ N, a best nonnegative rank-r approximation always exist in the sense that the following infimum inf{d (A, X) | X ∈ Ω 0 , rank + (X) ≤ r} is attained by some nonnegative tensor X r ∈ Ω 0 , rank + (X r ) ≤ r. In particular, the nonnegative tensor approximation problem X r ∈ argmin rank + (X)≤r A -X 1 Whenever possible, we will use up, vp, . . . , zp instead of the more cumbersome u

(1) p , u (2)
p , . . . , u (k) p to denote the vector factors in an outer product. It is to be understood that there are k vectors in "up, vp, . . . , zp," where k ≥ 3.

is well-posed. Here rank + (X) denotes the nonnegative rank of X and will be formally introduced in Section 4.

Tensors as hypermatrices

Let V 1 , . . . , V k be real vector spaces of dimensions d 1 , . . . , d k respectively. An element of the tensor product

V 1 ⊗ • • • ⊗ V k is called an order-k tensor. Up to a choice of bases on V 1 , . . . , V k , such a tensor may be represented by a d 1 × • • • × d k array of real numbers 2 , (3) A = a j 1 •••j k d 1 ,...,d k j 1 ,...,j k =1 ∈ R d 1 ו••×d k .
Gelfand, Kapranov, and Zelevinsky called such coordinate representations of abstract tensors hypermatrices [START_REF] Gelfand | Discriminants, resultants, and multidimensional determinants[END_REF]. It is worth pointing out that an array is just a data structure but like matrices, hypermatrices are more than mere arrays of numerical values. They are equipped with algebraic operations arising from the algebraic structure of

V 1 ⊗ • • • ⊗ V k :
• Addition and Scalar Multiplication:

For a j 1 •••j k , b j 1 •••j k ∈ R d 1 ו••×d k and λ, µ ∈ R, (4)
λ a j 1 •••j k + µ b j 1 •••j k = λa j 1 •••j k + µb j 1 •••j k ∈ R d 1 ו••×d k .
• Outer Product Decomposition:

Every A = a j 1 •••j k ∈ R d 1 ו••×d k may be decom- posed as (5) A = r p=1 λ p u p ⊗ v p ⊗ • • • ⊗ z p , a j 1 •••j k = r p=1 λ p u pj 1 v pj 2 • • • z pj k , with λ p ∈ R, u p = [u p1 , . . . , u pd 1 ] ⊤ ∈ R d 1 , . . . , z p = [z p1 , . . . , z pd k ] ⊤ ∈ R d k , p = 1, . . . , r.
The symbol ⊗ denotes the Segre outer product: For vectors x = [x 1 , . . . , x l ] ⊤ ∈ R l , y = [y 1 , . . . , y m ] ⊤ ∈ R m , z = [z 1 , . . . , z n ] ⊤ ∈ R n , the quantity x ⊗ y ⊗ z, is simply the 3hypermatrix x i y j z k l,m,n i,j,k=1 ∈ R l×m×n , with obvious generalization to an arbitrary number of vectors.

It follows from (4

) that R d 1 ו••×d k is a vector space of dimension d 1 • • • d k .
The existence of a decomposition (5) distinguishes R d 1 ו••×d k from being merely a vector space by endowing it with a tensor product structure. While as real vector spaces, R l×m×n (hypermatrices), R lm×n , R ln×m , R mn×l (matrices), and R lmn (vectors) are all isomorphic, the tensor product structure distinguishes them. Note that a different choice of bases on

V 1 , . . . , V k would lead to a different hypermatrix representation of elements in V 1 ⊗• • •⊗V k .
So strictly speaking, a tensor and a hypermatrix are different in the same way a linear operator and a matrix are different. Furthermore, just as a bilinear functional, a linear operator, and a dyad may all be represented by the same matrix, different types of tensors may be represented by the same hypermatrix if one disregards covariance and contravariance. Nonetheless the term 'tensor' has been widely used to mean a hypermatrix in the data analysis communities (including bioinformatics, computer vision, machine learning, neuroinformatics, pattern recognition, signal processing, technometrics), and we will refrain from being perverse and henceforth adopt this naming convention. For the more pedantic readers, it is understood that what we call a tensor in this article really means a hypermatrix.

A non-zero tensor that can be expressed as an outer product of vectors is called a rank-1 tensor. More generally, the rank of a tensor 2 The subscripts and superscripts will be dropped when the range of j1, . . . , j k is obvious or unimportant.

A = a j 1 •••j k d 1 ,...,d k j 1 ,...,j k =1 ∈ R d 1 ו••×d k , denoted
We use double brackets to delimit hypermatrices. rank(A), is defined as the minimum r for which A may be expressed as a sum of r rank-1 tensors [START_REF] Hitchcock | The expression of a tensor or a polyadic as a sum of products[END_REF][START_REF] Hitchcock | Multiple invariants and generalized rank of a p-way matrix or tensor[END_REF], [START_REF] Bro | Least squares algorithms under unimodality and non-negativity constraints[END_REF] rank(A) := min r A = r p=1

λ p u p ⊗ v p ⊗ • • • ⊗ z p .
The definition of rank in [START_REF] Bro | Least squares algorithms under unimodality and non-negativity constraints[END_REF] agrees with the definition of matrix rank when applied to an order-2 tensor.

The Frobenius norm or F -norm of a tensor

A = a j 1 •••j k d 1 ,...,d k j 1 ,...,j k =1 ∈ R d 1 ו••×d k is defined by (7) A F = d 1 ,...,d k j 1 ,...,j k =1 |a j 1 •••j k | 2 1 2 .
The F -norm is by far the most popular choice of norms for tensors in data analytic applications. However when A is nonnegative valued, then there is a more natural norm that allows us to interpret the normalized values of A as probability distribution values, as we will see in the next section. With this in mind, we define the E-norm and G-norm by ( 8)

A E = d 1 ,...,d k i 1 ,...,i k =1 |a j 1 •••j k | and A G = max{|a j 1 •••j k | | j 1 = 1, . . . , d 1 ; . . . ; j k = 1, . . . , d k }.
Observe that the E-, F -, and G-norms of a tensor A are simply the l 1 -, l 2 -, and l ∞ -norms of A regarded as a vector of dimension

d 1 • • • d k .
Furthermore they are multiplicative on rank-1 tensors in the following sense:

u ⊗ v ⊗ • • • ⊗ z E = u 1 v 1 • • • z 1 , (9)
u ⊗ v ⊗ • • • ⊗ z F = u 2 v 2 • • • z 2 , u ⊗ v ⊗ • • • ⊗ z G = u ∞ v ∞ • • • z ∞ .
The F -norm has the advantage of being induced by an inner product on

R d 1 ו••×d k , namely, (10) A, B = d 1 ,...,d k j 1 ,...,j k =1 a j 1 •••j k b j 1 •••j k .
As usual, it is straightforward to deduce a Cauchy-Schwarz inequality

| A, B | ≤ A F B F ,
and a Hölder inequality

| A, B | ≤ A E B G .
Many other norms may be defined on a space of tensors. For any 1 ≤ p ≤ ∞, one may define the l p -equivalent of ( 7), of which E-, F -, and G-norms are special cases. Another common class of tensor norms generalizes operator norms of matrices: For example if

A = a ijk ∈ R l×m×n and A(x, y, z) := l,m,n i,j,k=1
a ijk x i y j z j k denotes the associated trilinear functional, then A p,q,r := sup

x,y,z =0

|A(x, y, z)| x p y q z r defines a norm for any 1 ≤ p, q, r ≤ ∞. Nevertheless all these norms are equivalent (and thus induce the same topology) since the tensor product spaces here are finite-dimensional.

In particular, the results in this paper apply to any choice of norms since they pertain to the convergence of sequences of tensors.

The discussion in this section remains unchanged if R is replaced by C throughout (apart from a corresponding replacement of the Euclidean inner product in [START_REF] Comon | Independent component analysis: a new concept?[END_REF] by the Hermitian inner product) though a minor caveat is that the tensor rank as defined in [START_REF] Bro | Least squares algorithms under unimodality and non-negativity constraints[END_REF] depends on the choice of base fields (see [START_REF] Silva | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF] for a discussion).

Nonnegative decomposition of nonnegative tensors

We will see that a finite collection of discrete random variables satisfying both the naïve Bayes hypothesis and the Ockham principle of parsimony have a joint probability distribution that, when regarded as a nonnegative tensor on the probability simplex, decomposes in a nonnegative rank-revealing manner that parallels the matrix singular value decomposition. This generalizes Hofmann's probabilistic variant [START_REF] Hofmann | Probabilistic Latent Semantic Indexing[END_REF] of latent semantic indexing (lsi), a well-known technique in natural language processing and information retrieval that Harshman played a role in developing [START_REF] Deerwester | Indexing by latent semantic analysis[END_REF]. Nonnegative tensor decompositions were first studied in the context of parafac with nonnegativity constraints by the technometrics communities [START_REF] Bro | A fast non-negativity constrained least squares algorithm[END_REF][START_REF] Bro | Least squares algorithms under unimodality and non-negativity constraints[END_REF][START_REF] Carroll | Fitting of the latent class model via iteratively reweighted least squares candecomp with nonnegativity constraints[END_REF][START_REF] Krijnen | Contrastvrije oplossingen van het candecomp/parafacmodel[END_REF][START_REF] Paatero | A weighted non-negative least squares algorithm for three-way parafac factor analysis[END_REF]. The interpretation as a naïve Bayes decomposition of probability distributions into conditional distributions was due to Garcia, Stillman, and Sturmfels [START_REF] Garcia | Algebraic geometry of Bayesian networks[END_REF] and Sashua and Hazan [START_REF] Shashua | Non-negative tensor factorization with applications to statistics and computer vision[END_REF]. It is perhaps worth taking this opportunity to point out a minor detail that had somehow been neglected in [START_REF] Garcia | Algebraic geometry of Bayesian networks[END_REF][START_REF] Shashua | Non-negative tensor factorization with applications to statistics and computer vision[END_REF]: the naïve Bayes hypothesis is not sufficient to guarantee a nonnegative rank-revealing decomposition, one also needs the Ockham principle of parsimony, i.e. the hidden variable in question has to be minimally supported.

A tensor

A = a j 1 •••j k d 1 ,...,d k j 1 ,...,j k =1 ∈ R d 1 ו••×d k is nonnegative, denoted A ≥ 0, if all a j 1 •••j k ≥ 0. We will write R d 1 ו••×d k + := {A ∈ R d 1 ו••×d k | A ≥ 0}. For A ≥ 0, a nonnegative outer- product decomposition is one of the form (11) A = r p=1 δ p u p ⊗ v p ⊗ • • • ⊗ z p
where δ p ≥ 0 and u p , v p , . . . , z p ≥ 0 for p = 1, . . . , r. It is clear that such a decomposition exists for any A ≥ 0. The minimal r for which such a decomposition is possible will be called the nonnegative rank. For A ≥ 0, this is denoted and defined via rank

+ (A) := min r A = r p=1 δ p u p ⊗ v p ⊗ • • • ⊗ z p , δ p , u p , v p , . . . , z p ≥ 0 for all p .
Let ∆ d denote the unit d-simplex, i.e. the convex hull of the standard basis vectors in R d+1 . Explicitly,

∆ d := d+1 p=1 δ p e p ∈ R d+1 d+1 p=1 δ p = 1, δ 1 , . . . , δ d+1 ≥ 0 = {x ∈ R d+1 + | x 1 = 1}.
For nonnegative valued tensors, the E-norm has the advantage that (8) reduces to a simple sum of all entries. This simple observation leads to the following proposition stating that the decomposition in (11) may be realized over unit simplices if we normalize A by its E-norm.

Proposition 4.1. Let A ∈ R d 1 ו••×d k + be a nonnegative tensor with rank + (A) = r. Then there exist δ = [δ 1 , . . . , δ r ] ⊤ ∈ R r + , u p ∈ R d 1 -1 + , v p ∈ R d 2 -1 + , . . . , z p ∈ R d k -1 + , p = 1, . . . , r, where δ 1 = A E and u p 1 = v p 1 = • • • = z p 1 = 1, such that (12) A = r p=1 δ p u p ⊗ v p ⊗ • • • ⊗ z p .
Proof. If A = 0, this is obvious. So we will suppose that A = 0. By the minimality of r = rank + (A), we know that u p , v p , . . . , z p in ( 12) are all nonzero and we may assume that

u p 1 = v p 1 = • • • = z p 1 = 1
since otherwise we may normalize

ûp = u p / u p 1 , vp = v p / v p 1 , . . . , ẑp = z p / z p 1 ,
and set δp =

δ p u p 1 v p 1 • • • z p 1 ,
and still have an equation of the form in [START_REF] Silva | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF]. It remains to show that

δ 1 = A E .
Note that since all quantities involved are nonnegative,

A E = r p=1 δ p u p ⊗ v p ⊗ • • • ⊗ z p E = r p=1 δ p u p ⊗ v p ⊗ • • • ⊗ z p E .
By [START_REF] Carroll | Fitting of the latent class model via iteratively reweighted least squares candecomp with nonnegativity constraints[END_REF], the rhs can be further simplified to

r p=1 δ p u p 1 v p 1 • • • z p 1 = r p=1 δ p = δ 1 ,
as required.

Note that the conditions on the vectors imply that they lie in unit simplices of various dimensions:

(13) u 1 , . . . , u r ∈ ∆ d 1 -1 , v 1 , . . . , v r ∈ ∆ d 2 -1 , . . . , z 1 , . . . , z r ∈ ∆ d k -1 .
For k = 2, the above decomposition is best viewed as a parallel to the singular value decomposition of a matrix A ∈ R m×n , which is in particular an expression of the form

(14) A = r p=1 σ p u p ⊗ v p ,
where r = rank(A),

σ 2 = r p=1 |σ p | 2 1 2 = A F , and
u p 2 = v p 2 = 1,
for all p = 1, . . . , r. Here σ = [σ 1 , . . . , σ r ] ⊤ ∈ R r is the vector of nonzero singular values of A. If A is normalized to have unit F -norm, then all quantities in (14) may be viewed as living in unit spheres of various dimensions:

A ∈ S mn-1 , σ ∈ S r-1 , u 1 , . . . , u r ∈ S m-1 , v 1 , . . . , v r ∈ S n-1 where S d-1 = {x ∈ R d | x 2 = 1}
is the unit sphere in R d . For k = 2, the nonnegative matrix decomposition in Proposition 4.1 is one where the unit spheres are replaced by unit simplices and the l 2 -and F -norms replaced by the l 1 -and E-norms. An obvious departure from the case of svd is that the vectors in ( 13) are not orthogonal. Henceforth when we use the terms ntf and nmf, we will mean a decomposition of the type in Proposition 4.1. For a nonnegative tensor with unit

E-norm, A ∈ ∆ d 1 •••d k -1 , the decomposition in Proposition 4.1 has a probabilistic interpretation.
Let U, V, . . . , Z be discrete random variables and q(u, v, . . . , z) = Pr(U = u, V = v, . . . , Z = z) be their joint probability distribution. Suppose U, V, . . . , Z satisfy the naïve Bayes hypothesis, i.e. they are conditionally independent upon a single hidden random variable Θ. Let q 1 (u | θ), q 2 (v | θ), . . . , q k (z | θ) denote respectively the marginal probability distributions of U, V, . . . , Z conditional on the event Θ = θ. Then the probability distributions must satisfy the relation [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF] q(u, v, . . . , z) =

r θ=1 δ(θ) q 1 (u | θ)q 2 (v | θ) • • • q k (z | θ)
where δ(θ) = Pr(Θ = θ). Since the discrete random variables U, V, . . . , Z may take d 1 , d 2 , . . . , d k possible values respectively, the Bayes rule in [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF] can be rewritten as the tensor decomposition in [START_REF] Silva | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF], provided we 'store' the marginal distributions q 1 (u | θ), q 2 (v | θ), . . . , q k (z | θ) in the vectors u θ , v θ , . . . , z θ respectively. The requirement that r = rank + (A) corresponds to the Ockham principle of parsimony: that the model ( 15) be the simplest possible, i.e. the hidden variable Θ be minimally supported.

For the case k = 2, ( 15) is Hofmann's plsi [START_REF] Hofmann | Probabilistic Latent Semantic Indexing[END_REF], a probabilistic variant of latent semantic indexing [START_REF] Deerwester | Indexing by latent semantic analysis[END_REF]. While it is known [START_REF] Gaussier | Relation between plsa and nmf and implications[END_REF] that the multiplicative updating rule for nmf with kl divergence in [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF] is equivalent to the use of em algorithm for maximum likelihood estimation of plsi in [START_REF] Hofmann | Probabilistic Latent Semantic Indexing[END_REF], this is about the equivalence of two algorithms (em and multiplicative updating) applied to two approximation problems (maximum likelihood of plsi and minimum kl divergence of nmf). Since the em algorithm and the nmf multiplicative updating rules are first-order methods that can at best converge to a stationary point, saying that these two algorithms are equivalent for their respective approximation problems does not imply that the respective models are equivalent. The preceding paragraph states that the probabilistic relational models behind plsi and ntf (and therefore nmf) are one and the same -a collection of random variables satisfying the naïve Bayes assumption with respect to a parsimonious hidden variable. This is a statement independent of approximation or computation.

Nonexistence of globally optimal solution for real and complex tensor approximations

A major difficulty that one should be aware of is that the problem of finding a best rank-r approximation for tensors of order 3 or higher has no solution in general. There exists

A ∈ R d 1 ו••×d k such that (16) inf A - r p=1 λ p u p ⊗ v p ⊗ • • • ⊗ z p
is not attained by any choice of λ p , u p , v p , . . . , z p , p = 1, . . . , r. It is also in general not possible to determine a priori if a given A ∈ R d 1 ו••×d k will fail to have a best rank-r approximation. This problem is more widespread than one might imagine. It has been shown in [START_REF] Silva | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF] that examples of this failure happens over a wide range of dimensions, orders, ranks, and for any continuous measure of proximity (thus including all norms and Brègman divergence). Moreover such failures can occur with positive probability and in some cases with certainty, i.e. where the infimum in ( 16) is never attained. This phenomenon also extends to symmetric tensors [START_REF] Comon | Symmetric tensors and symmetric tensor rank[END_REF]. This poses some serious conceptual difficulties -if one cannot guarantee a solution a priori, then what is one trying to compute in instances where there are no solutions? We often get the answer "an approximate solution". But how could one approximate a solution that does not even exist in the first place? Conceptual issues aside, this also causes computational difficulties in practice. Forcing a solution in finite precision for a problem that does not have a solution is an ill-advised strategy since a well-posed problem near to an ill-posed one is, by definition, ill-conditioned and therefore hard to compute. This ill-conditioning manifests itself in iterative algorithms as summands that grew unbounded in magnitude but with the peculiar property that the sum remains bounded. This was first observed by Bini, Lotti, and Romani [START_REF] Bini | Approximate solutions for the bilinear form computational problem[END_REF] in the context of arbitrary precision approximations (where this phenomenon is desirable). Independently Harshman and his collaborators Kruskal and Lundy [START_REF] Kruskal | How 3-MFA data can cause degenerate parafac solutions, among other relationships[END_REF] also investigated this phenomenon, which they called parafac degeneracy, from the perspective of model fitting (where it is undesirable). In Theorem 6.1, we will prove the cheerful fact that one does not need to worry about parafac degeneracy when fitting a nonnegative parafac model.

The first published account of an explicitly constructed example of parafac degeneracy appeared in a study of the complexity of matrix multiplication by Bini, Capovani, Lotti, and Romani [START_REF] Bini | O(n 2.7799 ) complexity for n × n approximate matrix multiplication[END_REF]. However their discussion was for a context entirely different from data analysis/model fitting and was presented in notations somewhat unusual. Until today, many remain unconvinced that the construction in [START_REF] Bini | O(n 2.7799 ) complexity for n × n approximate matrix multiplication[END_REF] indeed provides an explicit example of parafac degeneracy and continue to credit the much later work of Paatero [START_REF] Paatero | Construction and analysis of degenerate parafac models[END_REF]. The truth is that such constructions are well-known in algebraic computational complexity; in addition to [START_REF] Bini | O(n 2.7799 ) complexity for n × n approximate matrix multiplication[END_REF], one may also find them in [START_REF] Bini | Approximate solutions for the bilinear form computational problem[END_REF][START_REF] Bürgisser | Algebraic complexity theory[END_REF][START_REF] Knuth | The art of computer programming, 2: seminumerical algorithms[END_REF], all predating [START_REF] Paatero | Construction and analysis of degenerate parafac models[END_REF]. As a small public service3 , we will translate the original construction of Bini, Capovani, Lotti, and Romani into notations more familiar to the technometrics communities.

In [START_REF] Bini | O(n 2.7799 ) complexity for n × n approximate matrix multiplication[END_REF], Bini, Capovani, Lotti, and Romani gave an algorithm that can approximate to arbitrary precision the product of two n × n matrices and requires only O(n 2.7799 ) scalar multiplications. The key to their construction is the following triplet of matrices which at first glance seem somewhat mysterious:

U =   1 0 1 0 1 0 0 0 ε ε 1 1 0 1 0   , V =     ε 0 0 -ε 0 0 -1 0 1 0 0 0 0 0 ε 1 -1 1 0 1     , W =     ε -1 ε -1 -ε -1 ε -1 0 0 0 0 1 0 0 0 -ε -1 0 ε -1 1 0 0 0 -1     .
We will show that these matrices may be used to construct a sequence of tensors exhibiting parafac degeneracy. We will assume that U has a 4th row of zeros and so U, V, W ∈ R 4×4 . As usual, u ij , v ij , w ij will denote the (i, j)th entry of the respective matrices. Let n ≥ 4 and x 1 , x 2 , x 3 , x 4 ∈ R n (or C n ) be linearly independent vectors. For ε > 0, define

A ε := 4 j=1 4 i=1 u ij x i ⊗ 4 i=1 v ij x i ⊗ 4 i=1 w ij x i .
Observe that

A ε = (x 1 + x 3 ) ⊗ (εx 1 + x 4 ) ⊗ (ε -1 x 1 + x 4 ) + x 3 ⊗ (-x 2 -x 4 ) ⊗ ε -1 x 1 + x 1 ⊗ x 4 ⊗ (-ε -1 x 1 -ε -1 x 3 ) + (εx 2 + x 3 ) ⊗ (-εx 1 + x 2 ) ⊗ (ε -1 x 1 + x 2 ) + (x 1 + εx 2 ) ⊗ (εx 3 + x 4 ) ⊗ (ε -1 x 3 -x 4 ).

It is straight forward to verify that

lim ε→0 A ε = A where A = x 1 ⊗ x 1 ⊗ x 1 + x 1 ⊗ x 3 ⊗ x 3 + x 2 ⊗ x 2 ⊗ x 1 + x 2 ⊗ x 4 ⊗ x 3 + x 3 ⊗ x 2 ⊗ x 2 + x 3 ⊗ x 4 ⊗ x 4 .
Note that the sequence A ε exhibits parafac degeneracy: as ε → 0, each of the summands becomes unbounded in magnitude but A ε remains bounded (and in fact converges to A).

Regardless of whether

x 1 , x 2 , x 3 , x 4 ∈ R n or C n , it is clear that for all ε > 0, rank(A ε ) ≤ 5.
Furthermore, one may show that rank(A) = 6 over C and therefore rank(A) ≥ 6 over R (cf. remarks at the end of Section 3). In either case, A is an instance in R n×n×n or C n×n×n where the approximation problem in (1) has no solution for r = 5 -since inf rank(X)≤5 A -X = 0 and rank(A) ≥ 6 together imply that argmin rank(X)≤5 A -X = ∅.

Hence the construction in [1] also yields an explicit example of a best rank-r approximation problem (over R and C) that has no solution.

Existence of globally optimal solution for nonnegative tensor approximations

As we have mentioned in Section 2, nonnegativity constraints are often natural in the use of parafac. Empirical evidence from Bro's chemometrics studies revealed that parafac degeneracy was never observed when fitting nonnegative-valued data with a nonnegative parafac model. This then led Harshman to conjecture that this is always the case. The text of his e-mail had been reproduced in [START_REF] Lim | Optimal solutions to nonnegative parafac /multilinear nmf always exist[END_REF].

The conjectured result involves demonstrating the existence of global minima over a non-compact feasible region and is thus not immediate. Nevertheless the proof is still straightforward by the following observation: If a continuous real-valued function has a non-empty compact sublevel set, then it has to attain its infimum -a consequence of the extreme value theorem. This is essentially what we will show in the following proof for the nonnegative parafac loss function (in fact, we will show that all sublevel sets of the function are compact). We will use the E-norm in our proof for simplicity, the result for other norms then follows from the equivalence of all norms on finite-dimensional spaces. Essentially the same proof, but in terms of the more customary F -norm, appeared in [START_REF] Lim | Optimal solutions to nonnegative parafac /multilinear nmf always exist[END_REF]. We will follow the notations in Section 4.

Theorem 6.1. Let A ∈ R d 1 ו••×d k be nonnegative. Then inf A - r p=1 δ p u p ⊗ v p ⊗ • • • ⊗ z p E δ ∈ R r + , u p ∈ ∆ d 1 -1 , . . . , z p ∈ ∆ d k -1 , p = 1, . . . , r is attained. Proof. Recall that R n + = {x ∈ R n | x ≥ 0} and ∆ n-1 = {x ∈ R n + | x 1 = 1}. We define the function f : R r × (R d 1 × • • • × R d k ) r → R by (17) f (T ) := A - r p=1 δ p u p ⊗ v p ⊗ • • • ⊗ z p E
where we let T = (δ 1 , . . . , δ r ; u 1 , v 1 , . . . , z 1 ; . . . ; u r , v r , . . . , z r ) denote the argument of f . Let D be the following subset of R r × (R

d 1 × • • • × R d k ) r = R r(1+d 1 +•••+d k ) , D := R r + × (∆ d 1 -1 × • • • × ∆ d k -1 ) r .
Note that D is closed but unbounded. Let the infimum in question be µ := inf{f (T ) | T ∈ D}. We will show that the sublevel set of f restricted to D,

E α = {T ∈ D | f (T ) ≤ α}
is compact for all α > µ and thus the infimum of f on D must be attained. The set

E α = D ∩ f -1 (-∞, α]
is closed since f is continuous (by the continuity of norm). It remains to show that E α is bounded. Suppose the contrary. Then there exists a sequence

(T n ) ∞ n=1 ⊂ D with T n 1 → ∞ but f (T n ) ≤ α for all n. Clearly, T n 1 → ∞ implies that δ (n)
q → ∞ for at least one q ∈ {1, . . . , r}. Note that

f (T ) ≥ A E - r p=1 δ p u p ⊗ v p ⊗ • • • ⊗ z p E .
is attained by X * ∈ T ⊂ S. The remaining implications (b) ⇒ (c) ⇒ (d) ⇒ (a) are obvious.

In the language of [START_REF] Bürgisser | Algebraic complexity theory[END_REF][START_REF] Strassen | Gaussian elimination is not optimal[END_REF], this says that 'nonnegative border rank' coincides with nonnegative rank. An immediate corollary is that the E-norm in Theorem 6.1 and Proposition 6.2 may be replaced by any other norm. In fact we will see later that we may replace norms with more general measures of proximity.

Corollary 6.3. Let A = a j 1 •••j k ∈ R d 1 ו••×d k be nonnegative and • : R d 1 ו••×d k → [0, ∞) be an arbitrary norm. Then inf A - r p=1 δ p u p ⊗ v p ⊗ • • • ⊗ z p δ ∈ R r + , u p ∈ ∆ d 1 -1 , . . . , z p ∈ ∆ d k -1 , p = 1, . . . , r is attained.
Proof. This simply follows from the fact that all norms on finite dimensional spaces are equivalent and so induce the same topology on

R d 1 ו••×d k +
. So Proposition 6.2 holds for any norms. In particular, the statement (b) in Proposition 6.2 for an arbitrary norm • is exactly the result desired here. Corollary 6.3 implies that the parafac degeneracy discussed in Section 5 does not happen for nonnegative approximations of nonnegative tensors. There is often a simplistic view of parafac degeneracy as being synonymous to 'between component cancellation' and thus cannot happen for nonnegative tensor approximation since it is 'purely additive with no cancellation between parts' [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF][START_REF] Shashua | Non-negative tensor factorization with applications to statistics and computer vision[END_REF]. While it provides an approximate intuitive picture, this point of view is flawed since parafac degeneracy is not the same as 'between component cancellation'. There is cancellation in nx ⊗ y ⊗ z -(n + 1 n )x ⊗ y ⊗ z but the sequence exhibits no parafac degeneracy. Conversely, the sequence of nonnegative tensors

A n = 0 1 1 1/n 1 1/n 1/n 1/n 2 ∈ R 2×2×2
may each be decomposed nonnegatively as [START_REF] Gelfand | Discriminants, resultants, and multidimensional determinants[END_REF] A

n = A + 1 n B + 1 n 2 C with A, B, C ∈ R 2×2×2 given by A = 0 1 1 0 1 0 0 0 , B = 0 0 0 1 1 0 1 0 , C = 0 0 0 0 0 0 0 1 ,
and each may in turn be decomposed into a sum of rank-1 terms. While there is no 'between component cancellation' among these rank-1 summands, it is known [START_REF] Silva | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF] that the convergence lim n→∞ A n = A exhibits parafac degeneracy over R 2×2×2 or C 2×2×2 , where there are decompositions of A n different from the one given in [START_REF] Gelfand | Discriminants, resultants, and multidimensional determinants[END_REF] exhibiting parafac degeneracy. That such decompositions cannot happen over R 2×2×2 + is precisely the statement of Proposition 6.2, which we proved by way of Theorem 6.1.

Brègman divergences

In many applications, a norm may not be the most suitable measure of proximity. Other measures based on entropy, margin, spectral separation, volume, etc, are often used as loss functions in matrix and tensor approximations. Such measures may not even be a metric, an example being the Brègman divergence [START_REF] Brègman | A relaxation method of finding a common point of convex sets and its application to the solution of problems in convex programming[END_REF][START_REF] Dhillon | Matrix nearness problems using Brègman divergences[END_REF][START_REF] Iusem | Bregman distance," and "Bregman function[END_REF], a class of proximity measures that often have information theoretic or probabilistic interpretations. In the definition below, ri(Ω) denotes the relative interior of Ω, i.e. the interior of Ω regarded as a subset of its affine hull; • is any arbitrary norm on R d 1 ו••×d k -again the choice of which is immaterial since all norms induce the same topology on Ω. Definition 7.1. Let ∅ = Ω ⊆ R d 1 ו••×d k be a closed convex set. Let ϕ : Ω → R be continuously differentiable on ri(Ω) and strictly convex and continuous on Ω. The function

D ϕ : Ω × ri(Ω) → R defined by D ϕ (A, B) = ϕ(A) -ϕ(B) -∇ϕ(B), A -B is a Brègman divergence if (i) For any fixed A ∈ Ω, the sublevel set L α (A) = {X ∈ ri(Ω) | D ϕ (A, X) ≤ α} is bounded for all α ∈ R. (ii) Let (X n ) ∞ n=1 ⊂ ri(Ω) and A ∈ Ω. If lim n→∞ A -X n = 0, then lim n→∞ D ϕ (A, X n ) = 0. (iii) Let (X n ) ∞ n=1 ⊂ ri(Ω), A ∈ Ω, and (A n ) ∞ n=1 ⊂ Ω. If lim n→∞ A -X n = 0, lim sup n→∞ A n < ∞, lim n→∞ D ϕ (A n , X n ) = 0, then lim n→∞ A n -X n = 0.
Note that D ϕ (A, B) ≥ 0 and that D ϕ (A, B) = 0 iff A = B by the strict convexity of ϕ. However D ϕ need not satisfy the triangle inequality nor must it be symmetric in its two arguments. So a Brègman divergence is not a metric in general.

Brègman divergences are particularly important in nonnegative matrix and tensor decompositions [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF][START_REF] Shashua | Non-negative tensor factorization with applications to statistics and computer vision[END_REF]. In fact, one of the main novelty of nmf as introduced by Lee and Seung [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF] over the earlier studies in technometrics [START_REF] Carroll | Fitting of the latent class model via iteratively reweighted least squares candecomp with nonnegativity constraints[END_REF][START_REF] Krijnen | Contrastvrije oplossingen van het candecomp/parafacmodel[END_REF][START_REF] Paatero | Positive matrix factorization: A non-negative factor model with optimal utilization of error estimates of data values[END_REF] is their use of the Kullback-Leibler divergence [START_REF] Kullback | On information and sufficiency[END_REF] as a proximity measure4 . The kl divergence is defined for nonnegative matrices in [START_REF] Lee | Learning the parts of objects by nonnegative matrix factorization[END_REF] but it is straightforward to extend the definition to nonnegative tensors. For A

∈ R d 1 ו••×d k + and B ∈ ri(R d 1 ו••×d k + ), this is D KL (A, B) = d 1 ,...,d k j 1 ,...,j k =1 a j 1 •••j k log a j 1 •••j k b j 1 •••j k -a j 1 •••j k + b j 1 •••j k ,
where 0 log 0 is taken to be 0, the limiting value. It comes from the following choice of ϕ,

ϕ KL (A) = d 1 ,...,d k j 1 ,...,j k =1 a j 1 •••j k log a j 1 •••j k .
We note that Kullback and Liebler's original definition [START_REF] Kullback | On information and sufficiency[END_REF] was in terms of probability distributions. The version that we introduced here is a slight generalization. When A and B are probability distributions as in Section 4, then A E = B E = 1 and our definition reduces to the original one in [START_REF] Kullback | On information and sufficiency[END_REF],

D KL (A, B) = d 1 ,...,d k j 1 ,...,j k =1 a j 1 •••j k log a j 1 •••j k b j 1 •••j k .
In this case D KL (A, B) may also be interpreted as the relative entropy of the respective distributions.

It is natural to ask if the following analogous nonnegative tensor approximation problem for Brègman divergence will always have a solution:

(19) X r ∈ argmin{D ϕ (A, X) | X ∈ ri(Ω), rank + (X) ≤ r}.
Clearly, the problem cannot be expected to have a solution in general since ri(Ω) is not closed. ) since if we set x n = y n = z n = [1, n -1 ] ⊤ , then as n → ∞, D KL (A, x n ⊗ y n ⊗ z n ) = 1 n 3 → 0. This is simply a consequence of the way a Brègman divergence is defined and has nothing to do with any peculiarities of tensor rank, unlike the example discussed in Section 5. This difficulty may be avoided by posing the problem for any closed (but not necessarily compact) subset of ri(Ω). Proof. Recall that S r :

= {X ∈ R d 1 ו••×d k + | rank + (X) ≤ r}.
The statement is trivial if rank + (A) ≤ r. So we will also assume that rank + (A) ≥ r + 1. Let µ be the infimum in [START_REF] Harshman | Data preprocessing and the extended parafac model[END_REF] and let α > µ. By (i) in Definition 7.1, the sublevel set L α (A) is bounded and so its subset K ∩ S r ∩ L α (A) = {X ∈ K ∩ S r | D ϕ (A, X) ≤ α} must also be bounded. Note that K ∩ S r is closed. Since ϕ is continuously differentiable on ri(Ω), the function X → D ϕ (A, X) is continuous and so K ∩ S r ∩ L α (A) is also closed. Hence D ϕ (A, X) must attain µ on the compact set K ∩ S r ∩ L α (A).

As one can see from the proof, Proposition 7.2 extends to any other measure of proximity d(A, X) where the function X → d(A, X) continuous and coercive. Of course this is just a restatement of the problem, the bulk of the work involved is usually to show that the proximity function in question has those required properties.

Aside: norm-regularized and orthogonal approximations

We have often been asked about norm-regularized and orthogonal approximations of tensors that are not necessarily nonnegative. These approximation problems are useful in practice [START_REF] Comon | Independent component analysis: a new concept?[END_REF][START_REF] Harshman | Data preprocessing and the extended parafac model[END_REF][START_REF] Paatero | Construction and analysis of degenerate parafac models[END_REF]. Nevertheless these always have optimal solutions for a much simpler reason -they are continuous optimization problems over compact feasible set, so the existence of a global minima is immediate from the extreme value theorem (note that this is not the case for nonnegative tensor approximation). In the following, we will let A ∈ R d 1 ו••×dr , not necessarily nonnegative.

Recall that O(n, r), the set of n × r matrices (r ≤ n) with orthonormal columns, is compact in R n×r . If we impose orthonormality constraints on the normalized loading factors in (2), i.

  For example let A = e ⊗ e ⊗ e ∈ R 2×2×2 + where e = [1, 0] ⊤ , then inf{D KL (A, X) | X ∈ ri(R 2×2×2 + ), rank + (X) ≤ 1} = 0 cannot be attained by any x ⊗ y ⊗ z ∈ ri(R 2×2×2 +

Proposition 7 . 2 .

 72 Let Ω be a closed convex subset of R d 1 ו••×d k + and A ∈ Ω. Let D ϕ : Ω × ri(Ω) → R be a Brègman divergence. Then

( 20 )

 20 inf{D ϕ (A, X) | X ∈ K, rank + (X) ≤ r}is attained for any closed subset K ⊆ ri(Ω).

  e. [u 1 , . . . , u r ] ∈ O(d 1 , r), . . . , [z 1 , . . . , z r ] ∈ O(d k , r), then it follows that |λ p | ≤ A 2 F for p = 1, . . . , r, i.e. λ ∈ [-A F , A F ] r ⊂ R r .Since the parafac objective is continuous and we are effectively minimizing over the compact feasible region[-A F , A F ] r × O(d 1 , r) × • • • × O(d k , r),this shows that orthogonal parafac always has a globally optimal solution.

And also to fulfill, belatedly, an overdued promise made to Richard when he was preparing his bibliography on parafac degeneracy.

This brought back memories of the many intense e-mail exchanges with Harshman, of which one was about the novelty of nmf. His fervently argued messages will be missed.

Acknowledgements

We thank reviewers for their many helpful comments.

Since all terms involved in the approximant are nonnegative, we have

where the last two equalities follow from [START_REF] Carroll | Fitting of the latent class model via iteratively reweighted least squares candecomp with nonnegativity constraints[END_REF] and

The proof essentially shows that the function f is coercive -a real-valued function f is said to be coercive for minimization if lim x →+∞ f (x) = +∞ [START_REF] Brinkhuis | Optimization: insights and applications[END_REF]. This is a standard condition often used to guarantee that a continuous function on a noncompact domain attains its global minimum and is equivalent to saying that f has bounded sublevel sets. A minor point to note is that had we instead optimized over a sum of rank-1 terms u ⊗ v ⊗ • • • ⊗ z, the proof would fail because the vectors u, v, . . . , z may be scaled by non-zero positive scalars that product to 1, i.e.

So for example (nx) ⊗ y ⊗ (z/n) can have diverging loading factors even while the outerproduct remains fixed. We avoided this by requiring that u, v, . . . , z be unit vectors and having a δ that records the magnitude.

The following proposition provides four useful characterizations of the statement that the function rank + : R d 1 ו••×d k + → R is upper semicontinuous. This is the nonnegative rank equivalent of a similar result in [START_REF] Silva | Tensor rank and the ill-posedness of the best low-rank approximation problem[END_REF].

can be approximated arbitrarily closely by nonnegative tensors of strictly lower nonnegative rank, i.e.

since any X ′ ∈ S r \T must have A -X ′ > A while we know that δ ≤ A . By the compactness of T , there exists X * ∈ T such that A -X * = δ. So the required infimum Next, the regularization proposed in [START_REF] Paatero | Construction and analysis of degenerate parafac models[END_REF] is to add to the parafac objective terms proportional to the 2-norm of each loading factor, i.e. [START_REF] Hitchcock | The expression of a tensor or a polyadic as a sum of products[END_REF] A -

2 ).

From constrained optimization theory, we know that, under some regularity conditions, minimizing a continuous function f (x 1 , . . . , x k ) under constraints x i 2 = r i , i = 1, . . . , k, is equivalent to minimizing the functional f (x 1 , . . . ,

2 for appropriate ρ 1 , . . . , ρ k ∈ R. In a finite dimensional space, the sphere of radius r i is compact, and so is the feasible set defined by x i 2 = r i , i = 1, . . . , k, and thus f must attain its extrema. In the same vein, ( 21) is equivalent to an equality constrained optimization problem and so norm-regularized parafac always has a globally optimal solution. This approach may also be applied to regularizations other than the one discussed here.


Huy Tran

Carmina Angelica Perez-Romero

Teresa Ferraro

Cecile Fradin

Nathalie Dostatni

Mathieu Coppey
email: mathieu.coppey@curie.fr

Aleksandra M Walczak
email: awalczak@lpt.ens.fr

Affiliations

LiveFly -a toolbox for the analysis of transcription dynamics in live Drosophila embryos

Keywords: Drosophila melanogaster, real-time monitoring, transcription dynamics, image analysis, early development

We present the LiveFly toolbox for quantitative analysis of transcription dynamics in live Drosophila embryos. The toolbox allows users to process two-color 3D confocal movies acquired using nucleilabeling and the fluorescent RNA-tagging system described in the previous chapter and export the nuclei's position as a function of time, their lineages and the intensity traces of the active loci. The toolbox, which is tailored for the context of Drosophila early development, is semi-automatic, and requires minimal user intervention. It also includes a tool to combine data from multiple movies and visualize several features of the intensity traces and the expression pattern.

Introduction

The adaptation of the RNA fluorescent tagging system to fly embryos allows us to capture the temporal aspects of transcription [START_REF] Garcia | Quantitative Imaging of Transcription in Living Drosophila Embryos Links Polymerase Activity to Patterning[END_REF][START_REF] Lucas | Live imaging of Bicoid-dependent transcription in Drosophila embryos[END_REF]. This system makes use of interactions between a specific Coat Protein (e.g. MCP) and an RNA stem loop generally derived from viral systems (e.g. MS2 stem loop) to monitor the dynamics of nascent RNA at the active transcription loci (see previous Chapter). These constructs, such as the MCP-GFP [START_REF] Garcia | Quantitative Imaging of Transcription in Living Drosophila Embryos Links Polymerase Activity to Patterning[END_REF][START_REF] Lucas | Live imaging of Bicoid-dependent transcription in Drosophila embryos[END_REF] and PP7 [START_REF] Fukaya | Enhancer Control of Transcriptional Bursting[END_REF] systems have been successfully implemented in the Drosophila embryo and used to reveal several important features of gene expression in the early phase of development [START_REF] Garcia | Quantitative Imaging of Transcription in Living Drosophila Embryos Links Polymerase Activity to Patterning[END_REF][START_REF] Lucas | Live imaging of Bicoid-dependent transcription in Drosophila embryos[END_REF][START_REF] Fukaya | Enhancer Control of Transcriptional Bursting[END_REF][START_REF] Bothma | Dynamic regulation of eve stripe 2 expression reveals transcriptional bursts in living Drosophila embryos[END_REF][START_REF] Desponds | Precision of Readout at the hunchback Gene: Analyzing Short Transcription Time Traces in Living Fly Embryos[END_REF][START_REF] Ferraro | Transcriptional Memory in the Drosophila Embryo[END_REF].

Obtaining biological information from the unprecedented real-time observations of transcription dynamics using the MCP-MS2 system requires processing of the movies obtained from microscopy.

Many available tools [START_REF] Santinha | iCellFusion: Tool for fusion and analysis of live-cell images from time lapsemultimodal mycroscopy[END_REF][START_REF] Young | Measuring single-cell gene expression dynamics in bacteria using fluorescence time-lapse microscopy[END_REF] were developed for bacterial colonies and are not well adapted for the analysis of fly images. They lack the flexibility of switching between a specialized automatic mode, which is essential for efficient processing, to manual intervention mode, which are often required to maximize the amount of samples acquired from each time-consuming experiment.

In this chapter, we present the LiveFly toolbox, which provides an interface for the analysis of the MCP-MS2 movies acquired in the early nuclear cycles, from nuclear cycle (nc) 10 to nc 14, using the confocal microscope. From the movie's two colored channels, one for nuclear labeling and one for observations of nascent mRNA dynamics, the interface exports the nuclei's position, lineage and the intensity of the active loci as a function of time. The processing is semi-automatic: it requires user input only at the beginning of the process, but allows users to intervene for manual corrections if needed. It also includes a Visualizer to combine the data from multiple movies and visualize several features of the intensity traces and the pattern. The toolbox has been used to study the hunchback patterning system [START_REF] Driever | A gradient of bicoid protein in Drosophila embryos[END_REF][START_REF] Porcher | The Bicoid morphogen system[END_REF], reducing processing time for each movie from days to hours.

The LiveFly toolbox has three modules as shown in Fig. 1, forming a pipeline for movie processing and analysis. We present the functions and guide the reader through each module in the following sections.

Materials

1 System requirements

1-MATLAB 2013 or above, with the Image Processing Toolbox 2-The LiveFly toolbox: available for download from: https://github.com/huytran216/LiveFly_toolbox (see Note 1) 3-MATLAB toolbox for Bioformat: version 5.2.1 or higher, available for download from: https://www.openmicroscopy.org/bio-formats/downloads/. The Bioformat toolbox (zipped) is to be extracted to the Tool folder of the LiveFly's main folder.

Livefly inputs

1-A time-lapse 3D movie captured with a confocal microscope. The compatible file formats are *.tif, *.czi or *.lsm commonly used in time-lapse microscopy. The movie has two channels: the first red channel visualizes the nuclei marker (either NUP or Histone) and the second green channel visualizes the MCP-GFP/PP7 protein (see Note 2). The movies usually capture a thin layer of the embryo so that there are no overlapping nuclei on the same z-stack. 2-Input parameters: time resolution and XYZ pixel resolution, embryo orientation and the coordinates of the imaging region relative to the embryo anterior/posterior poles.

3-(Optional) List of nuclei IDs, time frames to be ignored for spot detection (see Note 3 and section 3.3).

Methods

Nuclei Segmentation

The segmentation module works with the "nuclear" signal of His-RFP or NUP-RFP acquired using the red channel. It should be noted that the fly nuclear signal is distinct from that of bacterial colonies: nuclei are generally homogeneous in shape and size and are far apart from one another. However, nuclei division is very synchronous and rapid (it takes ~30 s for the whole nuclei population to fully divide).

Therefore, nuclear signal processing requires very robust nuclear tracking while the segmentation algorithm can be kept simple.

We opted to combine popular cell segmentation methods for static images [START_REF] Young | Measuring single-cell gene expression dynamics in bacteria using fluorescence time-lapse microscopy[END_REF][START_REF] Häkkinen | CellAging: a tool to study segregation and partitioning in division in cell lineages of Escherichia coli[END_REF][START_REF] Otsu | A threshold Selection Method from Gray-Level Histograms[END_REF] with frame-byframe image registration for precise and efficient nuclei tracking. For fast computation, we perform nuclei segmentation on 2D images created by a maximum projection of the image stacks onto the Z-axis (optical axis).

Maximum projection

This step creates a 2D movie from the histone/NUP channel of the original 3D movie for nuclei segmentation. The image in each frame of the 2D movie is the maximum projection of the corresponding 3D image. Added to this, between every two consecutive frames, an image registration, which records the displacement of image regions from one frame to another, is computed. This registration will be used to predict the movement of nuclei in the subsequent nuclei tracking and lineage construction. 2-Specify the 3D movie file name and path (Line 4-6), the output 2D movie name and the frame range of interest (Line 8-9) (see Note 4).

3-Run the script (Press F5)

This process will create two files in the input folder: a maximum projection file of the nuclei in TIF format (hereby referred as RED.tif), and a MATLAB data file containing the frame registration information (named RED_align.mat).

Nuclei Segmentation

The nuclei segmentation process is performed on the Segmentation platform. The user interface layout is shown in Fig. 2. At any time, you can load and save the segmentation progress. The segmentation in each frame is automatic based on the specified parameters, followed by manual correction if necessary.

Step by step (see Note 5):

1-In the LiveFly toolbox's main folder, run Cell segmentation/LiveFly_SEGMENTATION.m to start the nuclei segmentation user interface.

5-Manually increase or decrease the size of all masks by one pixel by pressing BW+/BWrespectively (see Note 7).

6-Update the number of pixels needed to be added to the mask's size in the Auto-size section (use negative numbers to decrease the mask size). This will enable automatic mask size adjustment following the segmentation process. The Auto-size parameter is generally constant throughout the processing of the whole movie. 2-Check the lineage tree by selecting Menu> Track> Plot Lineage Tree (as shown in Fig. 3) for abnormalities in the timing of nuclear divisions. You can identify the frames and nuclear ID with abnormalities by clicking on the lineage tree.

3-Most abnormalities are caused by imperfect image registration and nuclei masks being too small (there is no overlap between a nucleus and its intended ancestor). Examine the frames and nuclei with abnormalities by selecting Menu> Track> Check Track and perform manual corrections with the provided tool (add/remove/resize masks) on them if needed (see Note 11).

4-Once all the abnormalities are corrected, rerun the tracking process.

5-Reassign the nuclear cycle by selecting Menu> Track> Reassign nuclear cycle. You will be asked to specify the frame range for each nuclear cycle. The input will be used to determine which nuclear cycle each nucleus belongs to.

Spot detection

The spot detection module extracts the fluorescent spot location and intensity from the 3D green channel. Spots are first identified using the threshold method (see below) and then fitted with a 2D Gaussian landscape, from which the spot intensity is extracted. This module will also combine the spot information and the nuclei segmentation result to create the movie's complete output file, which is used in subsequent analyses. 2-a Run the script. An input dialog will appear asking you to either enter the debug mode to determine th or proceed with the spot detection. Select "Yes" to enter the debug mode.

Post-analysis correction

The exported output file is generally ready for subsequent analysis. However, due to many reasons, such as nuclei moving in and out of the captured region, it may be necessary to crop the movie frame, ignore specific nuclei, or a whole nuclear cycle.

1-Copy the file Data analysis/correction.m to the output file folder.

2-Specify the output file name and possible exceptions (e.g. bad nuclei, frames, nuclear cycles) in this file. See the comments in the correction file for details (see Note 15).

Pattern visualizer

The pattern visualizer is a graphical user interface that aids users to merge data from multiple fly movies of the same strain and genetic construct and visualize them, either individually or together. The interface is shown Fig. 5

Initial analysis

1-Click on Load movies to load data from all the movies. This step will read the movie output files, and extract information about each nucleus, including its position, nuclear cycle, lineage, and the intensity traces of active transcription loci. The exceptions specified in the correction file are taken into account. The program will extract several features of the traces, such as the total spot intensity (ΣI), the mean spot intensity (μI), time of first spot appearance (tinit), the total spot appearance duration (tactive), etc. See <Menu>Help>Show feature description> for the feature list and description.

2-Select a range of proper interphase duration for each nuclear cycle (see Note 17) by selecting

Refine interphase.

3-Extract the feature pattern by clicking Extract feature. This step will fit the trace feature values along the embryo anterior-posterior (AP) axis in each movie with a Hill function, from which the maximum intensity level, the pattern's border position (coordinate of the half-maximum expression point) along the AP axis and the Hill coefficient are extracted. It is possible to impose the same Hill coefficient or maximum intensity level on all movies by selecting the options in Hill fitting option. Note that only patterns in the movies that are selected in tick boxes will be fitted.

4-Following this step, all the nuclei information and the Hill fitting results can be exported in an excel file by selecting Menu>Export xls.It is recommended to save the dataset after this step to save time in future analysis.

Pattern Viewer

Once the movie data have been imported and processed, the trace feature patterns can be viewed in individual embryos (using Show single movie feature buttons) (Example in Fig. 6) or the movie patterns can be merged together (using Show all movie feature buttons) (Example in Fig. 7).

Edit the tick boxes list at the bottom of the interface to select the features of interest and determine how to visualize their patterns.

Notes

Note 1. The README.md file contains a link to a 3D sample movie and the microscopy configurations, which can be used to familiarize oneself with the whole process. Note 5. Once the basics of the segmentation process is understood, it is advisable to familiarize oneself with the hotkeys. The segmentation tool is designed so that most frames can be segmented with hotkeys alone. Actions involving pointers, such as mask adding/removing, are required infrequently during mitosis. Note 6. There are 4 view modes, accessible from the menu or the hotkey X.

1-Mode 1 shows the original nuclei image.

2-Mode 2 shows the original image with the overlays of the nuclei masks in red.

3-Mode 3 shows the original image with nuclei masks as wireframes. This is useful when visualizing densely populating masks (e.g. nuclear cycle 14) 4-Mode 4 shows in wireframe the nuclei mask plus the previous frame masks as red overlay.

Mode 4 is useful to recognize new cells that appear during segmentation. If a nucleus is detected in the current frame but has no overlap with any nuclei in the previous frame, its mask will appear as a blue wireframe instead of white. Note 7. The nuclei real radius is generally a few pixels larger than the detected size. It is advisable to adjust the mask size to the real size to minimize errors in the subsequent spot detection process. Note 8. View mode 4 is useful to check for the displacement of the nuclei between frames. While the image registration should work properly throughout the cycle interphase, it may fail during mitosis due to rapid mitotic waves. In this case, manual registration, which is done by a few clicks, is needed. Note 9. This step can be done with the hotkey Space. Generally, the automatic segmentation parameters require few adjustments during the whole process. The whole nuclear cycle interphase can be analyzed by using Space, inspection, Space, inspection etc. Note 10. Given the movement of nuclei during the development phase, it is recommended that the movie time resolution is kept below 20 s per frame if the users want to track nuclei across mitosis.

However, if users are only interested in specific nuclear cycles, the time resolution can be as slow as 40 s per frame. Note 11. Due to the fast movement of nuclei during mitosis, nuclei masks between frames may not overlap with one another leading to bad nuclei tracking. In this case, you should increase the mask size (using BW+) until they overlap. View Mode 4 is very useful for this correction. Note 12. The first spots that are visually detected are typically bright spots, the threshold values should therefore be slightly lower than half of their peak intensity. Note that setting the value of th(1) and th(2) too low will increase the rate of false positive detections (there is no spot in the nucleus but the program detects that there is), resulting in detected spots with near-zero intensity. This does not affect the spot intensity traces drastically but it may be computationally more expensive, due to the high number of spots detected. If the thresholds are set too low, a warning about too many (over 500) detected spots will appear. When it happens, interrupt the program and increase th(1) and th [START_REF] Lucas | Live imaging of Bicoid-dependent transcription in Drosophila embryos[END_REF] slightly.

Note 13. The output file, organized by columns (features) and rows (time frame, nuclei ID), can be opened with any sheet editing program. Note 14. Open the movie Check_image.tif to see a frame-by-frame comparison between the maximum projection image with the segmented nuclei and fluorescent spots. Note 15. The correction file, though not strictly required for the analysis, is strongly recommended. If no exceptions are declared, only the movie output file name should be specified.

In this correction file, one can specify:

1-The XY frame of interest (with xlim_left, xlim_right, ylim_up, ylim_down variables): this is useful to discard nuclei at the border of the capture region, which pop in and out frequently during the movie.

2-Frames to be discarded (with rm variable): spots detected in these frames are ignored (nuclei are still segmented and tracked). Specifying these frames is useful if the MCP-GFP forms aggregates right before mitosis, which are sometimes falsely detected as active transcription loci.

3-Whether the movie begins or ends in the middle of a nuclear interphase (with leftcensored, rightcensored variables). These incomplete nuclear cycles will be ignored.

Note 16. The data file is saved in the MATLAB *.mat format. This file can directly be loaded into the MATLAB workspace instead of using the Visualizer.

Note 17. This step will retain the nuclei that are detected within a certain time range. It is helpful to remove nuclei that move in and out of the captured image region.

1 -

 1 In the LiveFly toolbox's main folder, open Cell segmentation/LiveFly_CREATE_PROJECTION.m with the MATLAB editor.

2 -

 2 Load the 2D maximum projection movie via Menu> File> Load image… 3-Specify the average nuclear size and intensity threshold for nuclear segmentation. You will see the preview of the segmentation result in the Mask Panel. 4-Click on Segment for automatic segmentation of the current frame with the specified parameters. You can examine the segmentation result by cycling through the View Modes using View>Cycle View (X) (see Note 6). The View Mode 4 is recommended to check both the segmentation result and the overlap of nuclei masks between frames (which is critical in nuclei tracking).

7 - 8 -

 78 Select Menu> Edit> Realign for image re-registration if necessary. You can opt for automatic registration or manual registration using the MATLAB built-in tool (see Note 8). Manually correct the masks if needed. You can add new nuclei masks or remove existing masks by accessing Menu> Edit 9-Apply the automatic segmentation and mask size adjustment to the next frame by selecting Menu> Edit> Propagate Next/Prev (see Note 9). Continue this frame-by-frame automatic segmentation until parameter adjustments and manual corrections are required. 3.1.3 Nuclei tracking In the nuclei tracking process, the program uses the image registration to determine a nuclei's ancestor in the previous frame. Two nuclei in two consecutive frames are considered ancestor-descendant if their masks, realigned based on the image registration, overlap (see Note 10). If this nucleus has only one descendant in the next frame, this descendant is treated as the same nucleus. If it has more than one descendant in the next frame, the two most overlapping descendants are considered as the two daughters of that nucleus. 1-Press Menu> Track to begin nuclei tracking across frames.

1 -

 1 In LiveFly toolbox's main folder, open Spot detection> main_Spot_detection.m in the MATLAB editor. Set the raw 3D movie and 2D nuclei movie path (line 7-9). 2-Determine the threshold th=[th(1) th(2)] for spot detection. Here, th(1) is the minimum absolute spot intensity and th(2) is the minimum absolute spot intensity above the local background level (determined using a median filter). These two values are used for the initial identification of the fluorescent spots. Only the brightest spot(s) located in the nucleus space are considered active transcription loci. The threshold values are determined as follows:

2- b

 b Select the frame in which you know that there is a spot in the next input dialog.2-c Zoom into the spot location and determine the threshold values, based on the two panels (see Note 12) (Fig.4) using the built-in MATLAB Zoom In and Data Cursor. 2-d Quit the script by pressing Ctrl+C. 3-Set the value of th(1) and th(2) in the editor. Set the frame range of interest by modifying both it_start and it_end (line 15-16). Run the script. Skip the debug mode this time. 4-The complete movie output file, containing frame-by-frame nuclei locations, lineages and its spot intensity, are stored in the subfolder table_summary of the 3D movie's directory (see Note 13). check Spot_detection/Header.txt for the output file's header. A MATLAB file (in the format _config.m) specifying the movie's time and spatial resolutions, and the parameters for spot detection is also generated. 5-Check the spot detection results by opening the newly generated movie, located in the subfolder <output_images> folder (see Note 14).

  in with examples of data. To access the visualizer, run Data analysis/ LiveFly_VISUALIZER.m 3.4.1 Managing the dataset 1-Create a new dataset by selecting Menu>File>Create dataset. 2-Select Add… to add a new movie to the dataset. A dialog will appear inquiring for the movie correction file. 3-Edit the movie list with Remove/Move Up/Move Down buttons. 4-Select the nuclear cycle of interest in each movie by editing the tick boxes. 5-save/load the dataset any time using Menu>File (see Note16)

Note 2 .

 2 The color of the nuclei or MCP channel is not necessarily fixed. The nuclei channel can be specified in the variable channel.(line 8) Note 3. The information on exceptions (e.g. bad nuclei, frames, nuclear cycles) is only needed after analysis of the raw movie and can be updated anytime. Note 4. It is possible to analyze only a subsection of the movie by specifying the frame range (it_start, it_end variables (line 35-36)). In general, incomplete nuclear cycles or the second half of nuclear cycle 14 (after cellularization) are ignored.

Figure 1 .

 1 Figure 1. Modules of the LiveFly toolbox. (A) The Nuclei segmentation module extracts from the

Figure 2 .

 2 Figure 2. Nuclei segmentation user interface. The interface allows users to observe the

Figure 3 .

 3 Figure 3. Example of a nuclei lineage tree. The y-axis denotes the frame number. Each data

Figure 4 :

 4 Figure 4: Determining the threshold value for spot detection. Example of a zoomed-in spot in a

Figure 5 .

 5 Figure 5. Pattern Visualizer. The interface allows the user to manage movies acquired with

Figure 6 .

 6 Figure 6. Examples of individual movie transcription patterns in different nuclear cycles:

Figure 7 .

 7 Figure 7. Examples of merged movie transcription patterns in nuclear cycle 13: (A) Example

Figure 1 :

 1 Figure 1 :

Figure 2 :

 2 Figure 2 :

Figure 3 :

 3 Figure 3 :

Figure 4 :

 4 Figure 4 :

Figure 5 :

 5 Figure 5 :

Figure 6 :

 6 Figure 6 :

Figure 7 :

 7 Figure 7 :

Acknowledgements

This work was supported by a Marie Curie MCCIG grant No. 303561 (AMW), PSL IDEX REFLEX (ND, AMW, MC), ANR-11-LABX-0044 DEEP Labex (ND), ANR-11-BSV2-0024 Axomorph (ND and AMW) and PSL ANR-10-IDEX-0001-02. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.