Canonical tensor model through data analysis -- Dimensions, topologies, and geometries --

The canonical tensor model (CTM) is a tensor model in Hamilton formalism and is studied as a model for gravity in both classical and quantum frameworks. Its dynamical variables are a canonical conjugate pair of real symmetric three-index tensors, and a question in this model was how to extract spacetime pictures from the tensors. We give such an extraction procedure by using two techniques widely known in data analysis. One is the tensor-rank (or CP, etc.) decomposition, which is a certain generalization of the singular value decomposition of a matrix and decomposes a tensor into a number of vectors. By regarding the vectors as points forming a space, topological properties can be extracted by using the other data analysis technique called persistent homology, and geometries by virtual diffusion processes over points. Thus, time evolutions of the tensors in the CTM can be interpreted as topological and geometric evolutions of spaces. We have performed some initial investigations of the classical equation of motion of the CTM in terms of these techniques for a homogeneous fuzzy circle and homogeneous two- and three-dimensional fuzzy spheres as spaces, and have obtained agreement with the general relativistic system obtained previously in a formal continuum limit of the CTM. It is also demonstrated by some concrete examples that the procedure is general for any dimensions and topologies, showing the generality of the CTM.


Introduction
How to formulate a consistent theory for quantum gravity is one of the major problems in fundamental physics. While general relativity and quantum mechanics are believed to be the correct theories in their applied areas in physics, quantization of general relativity using standard (perturbative) quantum field theoretical method is hard due to non-renormalizable divergences from small scale quantum fluctuations [1]. 1 One very promising direction to solve the issue is to formulate spacetime and matter fields in a more fundamental way than continuous spacetime and point-like objects. In such attempts, spacetime is considered to be an emergent entity generated by the dynamics of more fundamental degrees of freedom.
Among the various approaches to quantum gravity in line with the thoughts above, tensor models are of much interest. They were originally proposed [4,5,6] as a generalization of random matrix models, which were successful for describing two-dimensional quantum gravity, with the hope of obtaining consistent theories for quantum gravity in dimensions higher than two. While the original models suffer from some difficulties in computability, 2 improved models called colored tensor models were introduced [9], that enabled various analytical computations in what is called 1/N expansions [10]. The results seem to show that the emergent spaces in the colored tensor models are like branched polymers [10,11], 2D quantum gravity, or mixtures [12,13], far from macroscopic spacetimes or our actual spacetime.
On the other hand, there is a model of quantum gravity with a causal structure, called causal dynamical triangulation, that successfully generates macroscopic spacetimes [14,15]. This is in contrast with the corresponding Euclidean model, called dynamical triangulation, which is not successful in this regard [16,17]. This fact suggests the importance of causal treatment in quantum gravity, and one of the authors of this paper proposed a new type of tensor model, which we call canonical tensor model (CTM) [18,19]. This is formulated as a first-class constraint system in Hamilton formalism with a canonical conjugate pair of real symmetric three-index tensors as its dynamical variables. 3 Its first-class constraint algebra closely resembles that of the ADM formalism of general relativity, and there indeed exists a formal continuum limit in which they agree [21]. There are also other remarkable connections to general relativity: The N = 1 CTM 4 agrees with the mini-superspace treatment of general relativity [22], and the classical CTM in the formal continuum limit agrees with a general relativistic coupled system of gravity, a scalar field, and higher spin fields in the Hamilton-Jacobi formalism [23].
The formal continuum limit above is obtained by a formal replacement of the discrete values of the indices of the tensors, a = 1, 2, . . . , N , to a continuous one, x ∈ R D . Therefore, this formal continuum limit is assuming a classical continuous spacetime from the beginning and does not tell anything about how such a space may emerge from the (quantum) dynamics 1 However, see for instance [2] on the recent developments in the asymptotic safety program [3]. 2 For the original models there are no so-called 1/N expansions, as do exist for the matrix models. Recently, introducing a traceless condition [7] or a pair of symmetric tensors [8] has been proposed as possible resolutions. 3 For a concise review of the CTM, see for instance the review section in [20]. 4 N denotes the dimension of the vector space associated to the tensor indices. In other words, each index takes values from {1, 2, . . . , N }. of the theory. A clue to the last question has been obtained in our previous paper [20]: The wave function of the quantum CTM has strong peaks at values of the tensors symmetric under Lie-groups. Since we know that various symmetries are associated to our spacetime, this result is encouraging. Then the next question which naturally arises is how to interpret such preferred values of the tensors as spacetimes.
The first step to answer this question would be to establish the correspondence between tensors and spacetimes. To this end, we introduce two well-known techniques in data analysis to the CTM, and formulate a systematic procedure to extract topological and geometric properties of spaces held by the tensors. The first technique is called tensor-rank decomposition (or CP, etc.) [24,25,26,27]. This is a certain generalization of the singular value decomposition of a matrix, and decomposes a tensor into a number of vectors. By regarding the vectors as points and their mutual inner products as quantities featuring distance relations among points, one can obtain a space with topological and geometric properties extracted from a tensor. Here, topological properties are extracted through the second technique from topological data analysis called persistent homology [28]. 5 Geometric structure is extracted by virtual diffusion processes over points which are also often used in data analysis [33,34].
After introducing some notions and ideas, we consider a homogeneous fuzzy circle and homogeneous fuzzy two-and three-dimensional spheres to demonstrate the method. We study the time evolution of the tensors corresponding to these fuzzy spaces under the classical equation of motion of the CTM and interpret them as the evolution of spacetime by the extraction procedure mentioned above. We compare the results with the classical equation of motion of the general relativistic system derived in a formal continuum limit of the CTM in a former paper [23] and find good agreement. This paper is organized as follows. In Section 2, we review some elementary facts about the tensor-rank decomposition, and interpret the vectors obtained from the decomposition as points. In Section 3, we give a systematic method of constructing real symmetric three-index tensors of fuzzy spaces corresponding to ordinary continuous spaces with any dimensions and topologies. In Section 4, we introduce the notion of neighborhoods in terms of mutual inner products among vectors representing points in the sense of Section 2. In Section 5, we review persistent homology, a technique from topological data analysis, and demonstrate how one can apply it to the fuzzy spaces. In Section 6, we point out that the derivative expansion previously performed in the formal continuum limit mentioned above can be represented in the form of a continuous tensor-rank decomposition. Here the vectors of the decomposition are expressed with the scalar and metric fields of the general relativistic system corresponding to the CTM. In Section 7, we present a method of obtaining the values of scalar and metric fields by virtual diffusion processes over continuously existing points, based on the expressions obtained in Section 6. In Section 8, the method developed for continuous cases in the preceding sections is generalized to discrete cases, namely for finite N , and a method to characterize the local distance structures in fuzzy spaces is presented. In Section 9, the classical equation of motion of the CTM is applied to the real symmetric three-index tensors describing fuzzy spaces, and time-evolution is roughly described as an increasing process of number of points and mutual distances among points. In Section 10, by applying the extraction procedure mentioned above, detailed analysis of the time-evolutions of homogeneous fuzzy S 1 , S 2 and S 3 is performed and good agreement is obtained with the general relativistic system corresponding to the CTM. In Section 11, we explicitly construct the real symmetric three-index tensors for fuzzy spaces with various dimensions and topologies to demonstrate the absence of limitations of the procedure, showing the generality of the CTM. The last section is devoted to a summary and future prospects. In Appendix A, we show the algorithm of the C++ program we made and used for the tensor-rank decomposition.

Tensor-rank decomposition and notion of point
In this section we introduce the tensor-rank decomposition, also often called CP-decomposition [24,25,26,27], and use it to interpret a tensor as a collection of points which form a space. Throughout this paper, unless otherwise stated, we consider tensors which are real, symmetric, and of three-way 6 : where σ denotes arbitrary permutations of a, b, c, and the indices run from 1 to N . This particular choice of tensors is considered because we are interested in applying methods developed here to the CTM, which has a similar setup. It is however straightforward to generalize the contents of this section to other types of tensors. We also assume that the underlying vector space admits an O(N ) symmetry, which is the natural symmetry for real inner product spaces and is the kinematical symmetry of the CTM.
One may define a point by the simplest possible tensor. In the present case of a real symmetric three-way tensor, the simplest possibility is given by Using the O(N ) symmetry in the underlying vector space, the general form for a single point is given by where v is an N -dimensional real vector. This implies that arbitrary single points are equivalent under the O(N ) symmetry up to the sizes. The tensor of the form (3) is also called a rank-one tensor.
A space may be described by a collection of such single points, leading to a tensor of the form, A tensor represented by a sum of R rank-one tensors is called a rank-R tensor, for the smallest possible R with a given P . Representing a given tensor in such a sum has various names such as tensor-rank decomposition, rank-one tensor decomposition, CP-decomposition, etc. [24,25,26,27], and essentially generalizes the single value decomposition for matrices. The decomposition always exists with a finite R for finite N .
An important fact about the decomposition of a tensor in our usage is that the set of vectors in the decomposition of a tensor has sorts of uniqueness [35,36], and therefore a space can be represented by points unambiguously 7 in our actual applications, unless the rank is taken to be unnecessarily too large in the approximate tensor-rank decomposition [37], which appears below. This is different from the matrix case, because the vectors in the single value decomposition of a matrix always have a large continuous ambiguity. For example, the expression, , without changing M ab . There are other differences and subtleties in the decomposition of a tensor in comparison with the matrix case. A tensor may have other tensor-rank decompositions with different R and v, though there are some proven cases with uniqueness (or a partial one). Here the least value of R is called the rank of the tensor. The rank of a tensor depends on the base field (namely, real or complex numbers for instance) and whether each rank-one tensor in the decomposition is restricted to be symmetric or not. Since each term in (4) is a symmetric real rank-one tensor, R should be more precisely referred to as real symmetric rank, and the decomposition (4) as symmetric tensor-rank decomposition over the reals. Unless otherwise stated, the tensor-rank decomposition in this paper is always assuming the form (4) with real values, and we simply ignore these specifications for brevity.
A typical rank is defined by any rank such that the set of tensors having the rank has positive measure in the whole space of the tensors. This means that a given tensor can be approximated as closely as one likes with a finite probability by a tensor with such a typical rank. It is known that there exists only a single typical rank for complex symmetric tensors with given w, N , where w denotes the amount of ways (the amount of indices) of a tensor. This rank is called generic symmetric rank, which we here denote by R g , and is explicitly given by with the following exceptions: R g (2, N ) = N , and R g is given by increasing the above formula by one for (w, N ) = (3,5), (4,3), (4,4), (4,5). Here · denotes the ceiling function. This statement is from the Alexander-Hirschowitz theorem [38] (See [39,35] for more details.). The number on the right-hand side of (5) is called expected rank, because it can be obtained by the simple number counting of the degrees of freedom.
In the real case, however, there exist a number of typical ranks for given w, N , the least value of which agrees with the generic rank of the complex case (See [40] for more details.). So, the space of real tensors is divided into a number of subregions each of which has a certain typical rank. The formula for the typical ranks is not known for general w, N except for some specific cases. For example, the typical ranks are 2 and 3 for (w, N ) = (3, 2) (See [41] for a table for three-way real tensors.).
The notion of typical rank implies that a given tensor can be approximated as closely as one likes by the form (4) with a typical rank. However, due to the lack of a general formula for typical rank and a practical systematic procedure, the tensor-rank decomposition is to optimize vectors v i to approximate a given tensor as much as possible with the form (4) with a value of R. So, practically, what we obtain is an approximate tensor-rank decomposition, rather than an exact (4), where the error ∆P abc should be made as small as possible. The error ∆P abc can be made (numerically) vanish if one takes R large enough, but R cannot be taken unconditionally large in practical computations. This is not only because the optimization process takes longer time for large R with larger degrees of freedom, but also because for larger R it becomes more difficult to avoid rough decompositions which contain mutual cancellations of the rank-one components (See [37,42] for more details). Therefore there exist various uncertainties in the decomposition. Is R taken large enough? Are the vectors optimized? How much of an error is reasonable to allow?
These uncertainties introduce uncertainties in results and are potentially very harmful when actually doing computations. In our applications, however, reasonable results are obtained by taking R reasonably large to make errors sufficiently small and repeating the optimization procedure several times to choose the best set of vectors. Here, for the optimization, we made a C++ program which implements the greedy algorithm described in [42] with an additional constraint. The program is described in some detail in Appendix A. It is worth noting that despite the possible numerical problems, the tensor-rank decomposition is well defined, so at least in principle we have a good notion of points corresponding to a tensor.

Real symmetric three-way tensors corresponding to fuzzy spaces
Real symmetric three-way tensors may be used to describe spaces through the algebra of functions acting on these spaces [43,44,45,46,47]. In this section we describe a systematic method to construct such tensors from their corresponding algebra. This method is particularly useful in constructing such tensors corresponding to homogeneous spaces invariant under Lie-group symmetries. A requirement for such tensors is that they should be invariant under the symmetric properties of the corresponding homogeneous spaces.
The rough idea of fuzzy spaces is to specify a space in terms of the algebra of functions on it rather than a coordinate system. This would be in accord with the fact that the relevant objects in physics are fields on a space rather than a space itself. Let us consider first an ordinary continuous space R D . In this case functions can be used to label points, because a single point, say ω 0 , can be identified by providing a localized function 8 f ω 0 (ω) = δ D (ω − ω 0 ). Therefore considering all the independent functions, which are f ω 0 with ω 0 ∈ R D , gives the whole space. The algebra of functions, reflects the pointwise structure of the continuous space. To get to more interesting cases one can modify this structure in various ways. Well known are the non-commutative spaces in which the function algebras are taken to be non-commutative (and usually associative) [48].
We modify the algebra in a different way, by picking up only functions corresponding to lower frequency modes than a cut-off and ignore all the other higher frequency modes. In this case, functions cannot represent single points ω 0 anymore, and the space necessarily becomes "fuzzy". This truncation gives a finite number of functions f a (ω) (a = 1, 2, . . . , N ) for a compact space M with a coordinate ω. The simplest way to obtain an algebra of functions is to truncate the products of functions by ignoring higher frequency modes. Such an algebra has the form f a (ω)f b (ω) = P ab c f c (ω) with structure coefficients P ab c taking the original values in the full algebra of the continuum case, but the summation over the modes c is truncated by c ≤ N . This procedure gives a commutative non-associative algebra. Now the structure coefficients can be extracted by considering This procedure naturally defines a three-way tensor corresponding to a space with fuzziness. 9 By considering real functions for all f a , (7) gives a real symmetric three-way tensor representing a fuzzy space. For a homogeneous space, P abc should be invariant under its symmetry, and therefore the function set must be taken so that it forms a certain representation of the symmetry and P abc is an invariant tensor.
A comment is in order. Comparing with the decomposition (4), one notices that (7) is nothing but a tensor-rank decomposition of P abc , where the index i is replaced by a continuum one ω. This does not mean that a fuzzy space defined by (7) requires an infinite R with a continuous index. In fact, for a finite N , the tensor-rank decomposition of P abc can always be performed with a finite R. Therefore a compact fuzzy space defined by (7) is always represented by a finite number of points. If one takes a limit back to the continuum (i.e. N → ∞), R should indeed become infinite, which is considered in a formal continuum limit of the CTM in Section 6.
Another comment concerns the ordinary continuum space. Let us take a basis of real functions by the delta functions mentioned above. Then (7) is given by P ω 0 ω 1 ω 2 = δ D (ω 0 − ω 1 )δ D (ω 0 −ω 2 ). Thus the ordinary continuum space is described by a continuous fully diagonal three-way tensor.
As a concrete example, let us consider a homogeneous fuzzy two-sphere. As described above, we take a cut-off L and take the spherical harmonics with angular momenta ≤ L as a function set. Then the real symmetric three-way tensor corresponding to a homogeneous fuzzy two-sphere is given by Here (l, m) takes l = 0, 1, . . . , L, m = −l, −l + 1, . . . , l, andỸ lm are the real functions defined byỸ where Y lm are the spherical harmonics and the star represents taking a complex conjugation. 10 The exponential damping factor is to make the high frequency cut-off smoother, which turns out to result in better behaved systems, as is explained in Section 4. In the section we apply the tensor-rank decomposition to the P in (8), and obtain the geometric picture in Figure 1, which clearly represents a spherical object.
The above procedure is general enough to construct various real symmetric three-way tensors corresponding to fuzzy spaces. Another simple example is a homogeneous fuzzy S 1 , which can be obtained by considering a real basis for functions on a circle. In Section 11, we explicitly construct real symmetric three-way tensors corresponding to spaces with various dimensions and topologies, and also non-orientability. There we find that the tensor-rank decomposition leads to topological and geometric interpretations in agreement with the corresponding continuum spaces. This demonstration proves the generality of our construction and that real symmetric three-way tensors can in principle represent any kinds of spaces. This last fact is particularly important for the generality of the CTM, in which the tensors are real symmetric three-way. This is in contrast with the other Euclidean tensor models [4,5,6,9], in which the number of ways (i.e., the amount of indices) of tensors is supposed to be in accord with the dimension of building simplicial spaces one considers.

Notion of neighborhoods in fuzzy spaces
In this section, we introduce the notion of neighborhoods around points in fuzzy spaces. Let us assume that a tensor-rank decomposition (4) is obtained for a given real symmetric three-way 10 More explicitly, we use a formula, tensor. A point i is represented by a vector v i , as discussed in Section 2. Here we introduce the notion quite naively from the inner product, but in Section 7 this form is justified in the CTM via a virtual diffusion process.
Let us define the neighborhood of a point i by the following set of points: 11 where the repeated index a is assumed to be summed over. Hereafter this standard convention is implicitly assumed for the indices originated from the tensor indices, e.g., a but not i, j in (10). The paramater c in (10) determines the size of the neighborhood: For larger c, the neighborhood becomes smaller, and vice versa.
As an example, let us consider a fuzzy two-sphere with L = 5 defined in Section 3. The dimension of the vector space of P is N = (L + 1) 2 = 36. Taking the rank to be R = 72, one can obtain a tensor-rank decomposition of P within a 2 percent error 12 . The left of Figure 1 is a histogram of the values of the inner products v i a v j a for i, j = 1, 2, . . . , R. The rightmost bins around 0.8 are composed of the self-inner products v i a v i a . The middle ones around 0.4 are composed of the inner products between the nearest neighbor points. Most of the inner products concentrate in a small region around the origin, which means that most of the points are not in their mutual neighborhoods for c > 0. The physical meaning of this concentration is that the fuzzy space respects locality, which is indeed what we hope for if we make N sufficiently large. As can be seen in Figure 2, this concentration around the origin becomes larger when the size of the fuzzy space is bigger, as more points are not in their mutual neighborhoods. This aspect can also be quickly understood by the fact that the probability for two independent N -dimensional vectors to have a relative angle θ is proportional to sin N −1 θ, which is the surface volume on a unit sphere at an angle θ from a vector. This phenomenon is called the concentration of measure in mathematical literature [49]. The right of Figure 1 shows the neighborhood relations among the fuzzy space points, which are connected if two points i and j satisfy v i a v j a > 0.2. Here the cutoff value is chosen so that the middle bunch of bins around 0.4 representing the nearest neighbor connections is well included. The figure clearly shows that the P defined in Section 3 represents a discrete analogue of a continuous two-sphere through our procedure. The topological aspect is discussed more precisely in terms of persistent homology in Section 5.
Let us comment on the importance of the damping factor e −l 2 /L 2 in (9), which smoothens the cut-off. Figure 3 shows the histogram of the inner products obtained from the fuzzy twosphere without the damping factor for L = 5. Comparing with the left of Figure 1, one can see that the peak around the origin is broadened into the negative values. This situation can While the former, corresponding to the sharp cut-off case, has a long-range oscillatory behavior with both positive and negative values, the latter, corresponding to the case with the damping factor, has the fast exponential damping behavior with positive values only. The former aspects lead to two major difficulties. As is discussed in Section 7, we use a virtual diffusion process to extract geometrical information. However, the diffusion equation with negative coefficients can have very unusual behavior. For example, diffusion is usually a process of easing concentrations, but with negative coefficients it can have diverging behavior, which contradicts the assumptions in Section 7. Another difficulty is that the broad distribution around the origin implies that there exist a substantial number of pairs of points having non-negligible inner products. This means that the fuzzy space has strong non-local features violating locality of an ordinary continuous space. These problems may also be solved by taking L sufficiently large, but for doing actual numerical computations we have to restrict ourselves to finite tensors and the damping factor is very useful.

Persistent homology
In Section 4, we have introduced the notion of neighborhoods. This characterizes local topological structure of fuzzy spaces. Global topological structure is also of much interest. In this section, we introduce the notion called persistent homology [28] as a method to extract the homological structure of the fuzzy spaces.
Let us first assume that a distance d(·, ·) between any pair of points on a fuzzy space (after a tensor-rank decomposition) is given. How to construct such a distance is discussed in due course. Let us denote the set of points which represent the fuzzy space by V .
Let us introduce a family of abstract simplicial complexes, parameterized by u, associated to a fuzzy space, which is called a Vietoris-Rips stream 13 and is denoted by VR(V, u). The complex, VR(V, u), is defined as follows: • The vertex set is given by V .
• For vertices i and j, the edge [ij] is included in VR(V, u) if and only if d(i, j) ≤ u.
• A higher dimensional simplex is included in VR(V, u) if and only if all of its edges are.
Since the Vietoris-Rips stream has the obvious property that VR(V, u) ⊂ VR(V, u ) for u ≤ u , it is called a filtered simplicial complex. A filtered simplicial complex has the following functorial property: For u ≤ u , the inclusion i : Given such a stream of simplicial complexes with the above functorial property, one can follow the creation and annihilation of the elements in the homology group of VR(V, u) while changing u. Here the filtration parameter u roughly corresponds to the resolution of distances. When u is smaller than any of the distances between points, all the points are independent; there is no non-trivial topological structure. When u is increased, points get connected to one another, and there appear edges and higher-dimensional simplices, leading to some non-trivial topological structure. An element of the k-homology group H k (VR(V, u)) corresponds to a k-cycle which is not the boundary of a k + 1-cycle, and the dimension of this group (the Betti number) corresponds roughly to the amount of holes with k-dimensional boundaries. When u is changed to become larger than the size of such a hole, the hole is filled by simplices and is not visible in the homology group. Therefore such a hole can be represented by an interval [u start , u end ), which represents its creation and annihilation and is called a Betti interval. For a given point set with mutual distances, there exist Betti intervals of various lengths. Each of them in principle is directly associated with the data, but the ones with long lengths are considered to be the intrinsic feature of a fuzzy space. On the other hand, the shorter ones are not stable against small perturbations, depending much on details, and are rather regarded as noises. This summarizes the idea of persistent homology, which extracts a topological structure from a discrete set of points with distances.
We now describe a simple way to construct a distance function d(·, ·), also called a metric, which we can use in the analysis of persistent homology below. We define for fixed c, where the neighborhood N c (i) is defined in (10). A path between i ∈ V and j ∈ V is defined as an ordered sequence of points p(i, j) = (p 1 , . . . , p n )| p 1 =i,pn=j , where the points in a pair (p k , p k+1 ) are always in each other's neighborhood. The length of a path is given by the sum of the distances of individual links, L p := n−1 k=1 d(p k , p k+1 ). The distance of two points is defined as the length of the shortest path between them, d(i, j) := min{L p |p(i, j)}. If two points i, j are not connected by a union of neighborhoods, we say d(i, j) = ∞. This distance function is relatively simple, whereas a more sophisticated notion of distance is introduced in Section 7. This simple distance function however is easy to calculate, and seems to be applicable to extract the intrinsic topological structure of a fuzzy space. This is because long-lived stable Betti intervals are not affected by detailed choices of distances, while noisy short-lived Betti-intervals may be changed. Figure 4 shows the Betti intervals for the fuzzy two-sphere with L = 5 and R = 72, using the distance function defined above. We used a program named "Ripser" 14 to compute the socalled persistence barcodes of the Vietoris-Rips stream and imported the data in Mathematica to construct the images. Here the graph of Betti k represents the Betti intervals for the homology group H k (VR(V, u)), which is the aforementioned persistence barcode. For 0 ≤ u < 1, u is smaller than any of the distances between the points by construction, and there are no edges. Since H 0 represents the homology class of connected components, the number of Betti 0 intervals equals that of points or the rank R in this region of u. For 1 ≤ u, the points are all connected to form one component such that there exists one long interval for Betti 0 until the maximum distance on a fuzzy space, in this case 7. For 1 ≤ u < 4, there exists one interval for Betti 2 , which represents the existence of a hole with a two-dimensional boundary. It vanishes at u = 4, when the hole is filled by simplices. Betti 1 is vanishing throughout the range of u. Thus, the long-life structure is observed to be dim(H 0 ) = dim(H 2 ) = 1, dim(H 1 ) = 0, topologically agreeing with an ordinary two-sphere.
Another application of persistent homology is to determine the topological dimension of a fuzzy space. If the analysis of persistent homology implies dim(H k ) = 0 with a certain k, one can know that the topological dimension of the fuzzy space should not be less than k. However, this method is not so useful, because, for example, a ball has a finite topological dimension but vanishing homologies except H 0 . A more useful way is to consider a reference point, say p, and a collection of points within a certain range of distance from it, d min ≤ d(·, p) ≤ d max , and study its persistent homology. One would expect that the collection of points forms a sphere of dimension being D − 1, where D is the dimension of the fuzzy space. This yields a local definition of the topological dimension around a reference point, and if this is the same for any choice of reference point except special points such as those on boundaries, the fuzzy space can be considered to have a well-defined topological dimension. In Figure 5, an illustrative example is shown for the same fuzzy two-sphere as the previous one.
6 Tensor-rank decomposition in a formal continuum limit In [23] the authors discussed the correspondence between the CTM and a general relativistic system by performing a derivative expansion of P in a formal continuum limit of the CTM. In the paper the authors considered derivatives up to the fourth order to analyze the equation of motion of the metric and a scalar field up to the second order of their derivatives, while additional higher spin fields must be taken into account in higher orders. In the current discussion, however, we are interested in the metric and the scalar field themselves with no necessity for their derivatives, and it is sufficient to consider a derivative expansion up to the second order.
In the formal continuum limit the indices of the tensor are assumed to become continuous coordinates in R D : Furthermore, a locality condition is imposed, which says P xyz = 0 only if x ∼ y ∼ z. This rough locality condition was mathematically translated to the tensor becoming a distribution and may be given by a derivative expansion: Distributions are defined by their action on test functions under integration. The authors showed that up to second order the expansion can be written as where the β(x) and β µν (x), which is symmetric, are the expansion coefficients, and f (x) is an arbitrary test function. The β and β µν in the expansion are fields on the space R D . These fields contain the degrees of freedom of the CTM in the formal continuum limit, which corresponds to a general relativistic system. The relation between the β fields and the scalar field φ and the metric field g µν of the relativistic system has been found by analyzing the equations of motion of the CTM and is given by [23] β where g = Det(g µν ). The test functions and β's are not usual scalar functions but have nonvanishing density weights, and they were fixed from their transformation properties under spatial diffeomorphisms which are part of the continuum limit of the SO(N ) symmetry transformation of the CTM. The weights are ] as in (16). In particular, these weights are taken so that the weights associated to each index of P xyz are [g 1/4 ] and the integral for an index contraction d D x P xab P xcd is invariant under diffeomorphisms.
Let us now consider a continuous analogue of the tensor-rank decomposition. For this we assume the form The integration form d D ωβ −2 (ω) is chosen such that the weight of the integration form is 0 so the w x (ω) are of weight 0 in ω. The w x (ω) still has a weight in x of [g 1/4 (x)], because each index of P xyz must have this weight as explained above. From (14) one can see that the w x (ω) can also be given by a derivative expansion of the form where h.o. means higher orders, δ Here, as explained above, w x (ω) must have the weight of [g 1/4 (x)]. Let us use the convention that the density weight of the delta function, say δ D (x − y) having a total weight of [g 1/2 ], is equally distributed over both arguments x, y. Then the weights of the fields in (18)  By putting (18) into (17), multiplying test functions, integrating over their arguments, and comparing the result to (15), one finds According to the interpretation given in Section 2, the vector w x (ω) represents a single point labeled by ω. With these w x (ω) one can define a quantity similar to the Euclidean inner product between two points. Since the weight of w x (ω) is [g 1/4 (x)] an invariant quantity can be obtained by where Here we have used (16) to obtain the last field theoretical expressions, and it is apparent that the density weights provided by g −1/4 cancel the weights from the delta functions to make K(ω,ω) a scalar in ω andω.
As above, the quantity K(ω,ω) transforms as a scalar under diffeomorphisms on ω,ω. However, this invariant feature should be particular to the continuum case, because the (almost) uniqueness of the tensor-rank decomposition for finite N , mentioned in Section 2, does not allow such degeneracies of expressions in the discrete case. On the other hand, it would be important to use a corresponding similar form even in the discrete case, because it can be expected to converge to this physically meaningful invariant form in a continuum limit with N → ∞. Therefore we follow similar steps taking care of invariant forms as in this section, when we discuss a discrete analogue in Section 8.
The expression (21) of K(ω,ω) is given by an expansion in terms of the derivatives of delta functions, inheriting the locality imposed for P xyz below (13). The physical meaning of this fact is that locality is respected by the mutual relations among points. Similar inner products, v i a v j a , can be considered for the discrete case and characterize the local distance structures of a fuzzy space. This aspect of the inner products has already been used to define local neighborhoods around points in Section 4. In the following sections, this aspect is further pursued in more detail.

Distances by a virtual diffusion process
The main purpose of the present and the following sections is to find a notion of distances on fuzzy spaces in terms of a diffusion process by using the knowledge of the preceding section.
For this purpose we need to relate the distances given by the metric field in the continuous theory to the discrete case. This is done through K(ω,ω) defined by the inner product in (20), which can also be interpreted as a second order differential operator as shown in (21). This operator defines a virtual diffusion process on a continuum space, which can be easily replicated on a discrete space to extract corresponding continuum quantities.
Using virtual diffusion processes to interpret discrete systems similarly to continuous ones is common in the literature of data analysis and quantum gravity as they can be interpreted in a similar way for both the discrete and continuous cases. For instance in data analysis diffusion processes are often used in order to define distance functions in data sets [33,34]. In quantum gravity the use of diffusion processes is also appreciated [50], as they can be defined similarly for continuous and discrete spaces and allow one to construct well-defined observables such as the spectral dimension [51]. Our strategy is to consider a diffusion process defined by K(ω,ω) in (21) and extract the geometric and scalar field data from it. As defined in (20), K(ω,ω) is the continuous analogue of the inner product v i a v j a , so we can easily relate it to the discrete model and find a notion of distances there, which is done in the following section.
The diffusion equation we consider in the continuum case is given by where ρ is a scalar field representing the density of a virtual diffusing material, K(ω,ω) is given in (21), and β(ω) −2 makes the volume element to have weight zero. By using (21) and (22) and performing partial integrations, we obtain where In general the diffusion equation (24) cannot be solved analytically and one would have to rely on numerics. Rather than doing so, let us restrict ourselves to extracting only short distances by the diffusion equation. For this purpose, it is enough to consider a localized initial condition like ρ(ω, 0) = δ D (ω − ω 0 ) for arbitrary location ω 0 and a short time period of evolution 0 ≤ s 1. Since ρ is non-vanishing only in a small distance region around ω 0 under such a short period of time, one can regard B's as constants, assuming B's are smooth enough in ω. Then we can solve (24) and obtain where ρ 0 is an overall constant factor, and δω = ω − ω 0 . The expressions are still complicated to actually work with, but we can further assume the covariant derivatives of β's to vanish in the homogenous case and obtain Thus the diffusion process can determine a conformally rescaled metric e 2φ(ω) g µν (ω). 15 A comment is in order. We could have used d Dω g(ω) as the volume element in (23). In this case, (27) is changed to So in this case the φ and g µν fields completely decouple in the diffusion process. Namely, this choice gives more direct meaning to the coefficients of the diffusion equation from the point of view of the identification of the fields obtained in [23]. However, this choice is practically difficult to implement for the discrete case, since there is no natural way to know √ g while defining the kernel. Moreover, in the presence of a scalar field, there is no canonical way to take a particular choice of the metric from the ambiguity of the conformal transformation with the scalar field. This is what is called frame dependence, and various choices are possible depending on usages such as the Einstein frame normalizing the Einstein term. Therefore, we rather use β(ω) −2 for the volume element as above, which is much easier to implement in the discrete case, as is done in Section 8. Though the fields are not separated in the coefficients in this method, it gives a way to extract β and β µν /β in a straightforward manner and can equivalently determine φ and g µν through the relation (16).

Distances on fuzzy spaces
In this section, we discuss the actual process of determining distances between points on fuzzy spaces by considering discrete analogue of the method developed in Section 7. In fact, due to the difference between continuum and discrete spaces, we find an issue that there exist some offsets in the distances determined by the procedure for the discrete case. We propose a provisional solution to this issue, and get acceptable results in the actual application in Section 10. However, a more satisfactory resolution is desirable.
Let us start with the continuum case. One can determine distances between nearby 16 points by measuring the third term in (26). This third term is a damping function in s and generally tiny compared to the first term linear in s. Therefore it is practically (or numerically) difficult to measure, if the first term exists. To circumvent the situation, it is more convenient to replace the kernel in (23) by Then, assuming the homogenous case (27), the problematic first term disappears from (26), as well as the second term. Thus, we obtaiñ whereρ is the density function after the replacement (29). The maximum ofρ is located at s = s max satisfying Thus, by measuring s max of diffusion processes, one can determine the distance squares between arbitrary nearby points. Here note that the distances are defined with respect to the conformally rescaled metric e 2φ(ω) g µν (ω), as noted below (27). Now let us discuss the discrete case. The discrete analogue of the tensor-rank decomposition to the continuum (17) is given by and the diffusion kernel corresponding to (20) is given by A non-trivial part is how to determine β(i) from the tensor-rank decomposition. We consider a self-consistency condition given by This is derived from the following continuum counterpart, d Dω β(ω) −2 K(ω,ω) = 1, which can be proven for the homogeneous case due to the vanishing of the derivatives of β's. Since the relation between the two tensor-rank decompositions, (4) and (32), is given by the condition (34) can be rewritten as This condition determines β(i) from the decomposition (6), and hence w i a by (35). This process is used in our numerical analysis.
Next, let us consider the discrete version of the diffusion equation (23) with the substitution (29). This is given byK Here note that the last term in the first line is actually δ ij by using (34). Note also that K(i, j) has been defined so that the weight associated to i, j is 1/2 and the same forρ. This is to makeK(i, j) symmetric to simplify the following discussions. The other assignments of weights would be possible, like considering the diffusion equation with simpler assignments of weights, dρ(i)/ds = jK (i, j)β(j) −2ρ (j), whereρ andK have no weights. This weight assignment was actually used in the continuum discussions. However, this requires us to treat K(i, j)β(j) −2 , which is asymmetric and makes things non-obvious about eigenvalue problems and the symmetry of distances under mutual permutations of points. Therefore we employ the above symmetric assignment, which is indeed equivalent to any asymmetric assignment, because they are related by a similarity transformation β(i) wK (i, j)β(j) −w .
The distance square between arbitrary nearby points, i and j, can be determined from s max (i, j) at whichρ(i, s; j) takes the maximum value in s. BecauseK is symmetric,ρ(i, s; j) andρ(j, s; i) give the same distance, which guarantees the symmetry of the distances. 17 To study time evolutions of fuzzy spaces, we are interested in the time dependence of their sizes. For a given fuzzy space, one can perform the tensor-rank decomposition and solve the diffusion equation (37) to obtain s max (i, j) for arbitrary nearby points. This method cannot directly be used for long distances, because of the assumptions made in the derivation in Section 7. One would also think that the distance between two arbitrary points could be determined, even if it is large, by considering the shortest path connecting the two points, where its length is the sum of short lengths along the path. However, this turns out not to be justified because of the existence of the offsets explained below: The contribution of the offsets becomes considerable by being multiplied by the number of the short length portions along the path. Instead, since we are only considering homogeneous fuzzy spaces in this paper, we characterize the local distance structures of the fuzzy spaces and assume them to be proportional to their whole sizes.
To do so, let us first recall the distance introduced in Section 5, where we have discussed topological structures of fuzzy spaces. The distance is defined to take d(i, j) = 1, if two points i and j are in their local neighborhoods. Let us call this a topological distance and denote it by d t (i, j), because this is determined by topological relations of neighborhoods. Here, in the example of the fuzzy space in Figure 1, the two points i, j with d t (i, j) = 1 are those which have the inner products v i a v j a in the range between 0.4 and 0.6. Then topological distances d t (i, j) between any points i, j are defined by taking the shortest paths as given below (12). Now characteristic local distances of a homogeneous fuzzy space can be obtained by considering s max (r), which is an average value of s max (i, j) over all the pairs of i, j with d t (i, j) = r. This is plotted against r for the examples of fuzzy S 1 and S 2 in Figure 6. The data points in the figure are fitted with a quadratic function a 0 + a 1 r 2 with coefficients a 0 , a 1 . In the continuum case, a 0 = 0 and the distance is strictly proportional to √ s max , but in the present case, the offset is non-vanishing, a 0 = 0. This would be understandable because the diffusion process is from points to points at small s and it becomes continuous only after larger s. Therefore the discreteness is apparent in small s, and may generate such a difference from the continuous case. To regard a 0 negligible, we have to consider larger r such that s max (r) a 0 . On the other hand, we cannot take r too large because this violates the assumptions made in Section 7. In Section 10, we take r = 4 and obtain some acceptable results.

Time evolutions of fuzzy spaces in the CTM
As is explained, the equation of motion of the CTM gives a first-order differential equation in time for a real symmetric three-way tensor. In this section, we regard the solutions to the differential equation as time-evolutions of the tensors corresponding to fuzzy spaces, and study some of their elementary properties. One observation is that the time-evolutions increase the number of points forming fuzzy spaces starting from one, in the sense which will be described more precisely.
The degrees of freedom of the CTM are a pair of canonically conjugate real symmetric three-way tensors, which satisfy the fundamental Poisson brackets, where the summation is over all the possible permutations of d, e, f for the consistency with the permutation symmetry of the tensors. The classical equation of motion of the CTM is given by where X is Q or P , and the Hamiltonian H is given by a linear combination of the first-class constraints, H a and J ab , as Here N a and N ab are freely choosable generally time-dependent variables corresponding to the lapse function and the shift vector in the ADM formalism of general relativity. The explicit expressions of the constraints are given by In this paper, we put N ab = 0, since the corresponding term in (40) is just a generator of timedependent SO(N ) transformations, which are irrelevant if we are only interested in SO(N ) invariant quantities like the inner products v i a v j a . Then the equation of motion of P abc is given by dP abc dt = −N d (P dae P ebc + P dbe P eca + P dce P eab ) .
In this paper, we do not consider the equation of motion of Q abc , because the interpretation of the equation of motion in the continuum language (namely, general relativity) is only known for P abc [23]. We also do not consider a term, λ Q abb , which can be added to H a , because it causes an issue concerning locality in the classical equation of motion of P abc [21].
Let us consider the time evolution of the homogeneous fuzzy two-sphere defined in Section 3. The index set is given by a = (l, m), where l and m are integers satisfying 0 ≤ l ≤ L, −l ≤ m ≤ l with a cut-off L. Let us use 0 to represent the index (0, 0) for notational simplicity. We take the fuzzy two-sphere in (8) as the input of P abc at t = 0. For the lapse function, we take N 0 = 1 3 for convenience, where we also have to take N a = 0 for a = 0 to keep the SO(3) symmetry of the homogeneous fuzzy two-sphere under the time-evolution. 18 First of all, let us point out an important property general for homogenous cases. This has a similarity to the property first pointed out for N = 1 [22]. Due to the SO(3) symmetry, P 00a = 0 for a = 0. Then, by putting a = b = c = 0 in (42), one can find that the equation of motion of P 000 in (42) decouples from the others, and obtain dP 000 dt = −P 2 000 .
The solution is where P 000 (0) = 1 has been assumed as the normalization of the initial condition. As can be seen in (44), the solution diverges at t = −1 and monotonically decreases as t increases. The time t = −1 can be considered to be the time of birth of the fuzzy two-sphere as is explained below. The time-dependence of the other components can be computed numerically, and we used the Runge-Kutta method for the purpose.
Once a solution is obtained, one can perform the tensor-rank decomposition of P abc (t) at each t by the program described in Appendix A. Figure 7 plots the error ratio ∆P 2 /P 2 of the decomposition, where ∆P is the error of the approximate tensor-rank decomposition in (6). In Figure 7, P abc (0) is taken to be the fuzzy two-sphere in (8) with L = 5. The plot shows that the rank R must be increased for larger t, if one wants to keep the error ratio being suppressed under a certain value. In the t → −1 limit, one can numerically find that P 000 dominates over all the other components, meaning that P abc (t) approaches a rank-one tensor. This explains the rapid decaying behavior of the error ratio in the small t region. By regarding the rank to be equivalent to the number of points forming a space, the time evolution of P abc (t) can be regarded as that starting from one point at t = −1 and gradually increasing the number of points. 19 18 Physically, this is to consider the lapse which is uniform on the fuzzy two-sphere, or taking the space-like slices in which the evolution is described in a spatially uniform manner. An important thing here is that, because of the first-class nature of the constraints, this is just a gauge choice, but not a particular choice of a time-evolution. 19 This should not be considered as a mathematically rigid statement and is merely an approximate one, which would be relevant in practical physical applications containing errors or quantum fluctuations. One can easily see that, if an exact decomposition of P abc (t) is given at one time, the differential equation (42) can be rewritten in a closed form with the vectors v i a (t) only. This means that the rank of P abc (t) is constant in the course of changing t. Changing rank is needed only if a tensor-rank decomposition is approximate with an error and one wants to keep the error ratio under a certain value throughout a time-evolution. A more mathematically rigid statement is left for future study.  Figure 7: The error ratio ∆P 2 /P 2 of the tensor-rank decomposition of P abc (t) for a fuzzy two-sphere with L = 5 taken for P abc (0). P abc (t) with larger t requires higher ranks to keep the preciseness of the approximate tensor-rank decomposition. For R = 82, the error is numerically consistent with zero, suggesting that the actual rank of the tensor is 82. Another aspect related to the time evolution appears in the inner products v i a v j a . Figure 8 shows the histograms of the inner products v i a v j a at t = −0.75 and t = 25 for a homogeneous fuzzy two-sphere with N = 36 and R = 78. The inner products are shifted to the positive values for smaller t, meaning that the distances between points become shorter, and vice versa. This is consistent with the rough picture that the two-sphere becomes larger in time, as discussed above in the sense of the number of points. In the following section, we perform more detailed analysis with comparison with the general relativistic system derived in [23].

Correspondence to a general relativistic system
A discrete theory of quantum gravity is expected to reproduce a general relativistic system in a certain classical continuum limit. In [23], a general relativistic system corresponding to the CTM has been obtained by taking a formal continuum limit, which is to formally replace the discrete values taken by the indices, a = 1, 2, . . . , N , to the continuum coordinates, x ∈ R D . This formal replacement is obviously unsatisfactory from the view point of quantum gravity, because no considerations of dynamics to lead to the replacement were made. There is not only this difficult problem of dynamics of emergent macroscopic spaces, but also another related issue in this formal continuum limit that it is not given as any limit of N → ∞ but is rather given as a sudden formal replacement from the discrete to the continuum indices. The question we consider in this section is whether fuzzy spaces with large N can well be described by the continuum general relativistic theory obtained previously in [23] or not. We perform some detailed studies of the time evolutions of the homogeneous fuzzy S 1 , S 2 and S 3 by using the methods developed in the former sections and compare the results with the equation of motion of the general relativistic system. We obtain good agreement with the continuum theory at least for these homogeneous cases.
Let us first discuss the continuum side. The equations of motion of β(t, x) and β µν (t, x) in the continuum general relativistic theory are given by [23] 20 where n(t, x) is the lapse function, R µν is the Ricci curvature of the space, and we have ignored a number of covariant derivative terms in the original equations, since we are considering homogeneous spaces. In the case of homogenous spaces with uniform time evolutions n(t, x) = n(t), one can assume that β(t, x), β µν (t, x) are given by products of functions separately depending on time or space. Namely, we can write β µν (t, x) = β 2 (t)β µν (x) and β(t, x) = β(t), and the second equation in (45) can be rewritten as where c 1 is a constant proportional to the curvature on the space. The equation for β remains the same as in (45). Then, by taking n = 1/9 for convenience, one can obtain the solution to the above equations as where t 0 is supposed to be the time of birth of the space, β 0 2 is an integration constant, and c R = c 1 β 0 2 /6. As explained in Section 9, t 0 = −1 in our case. The solution leads to the following time dependence of B −1 µν by taking the inverse of the expression in (27): 20 The equations of motion are taken from Section VII of [23] with the consideration of the gauge condition ββ µν = g µν / √ g and the change of the allover minus sign for a convention. Let us compare the solution with the time evolutions of the fuzzy spaces. We consider homogeneous fuzzy S 1 with N = 31, R = 46, S 2 with N = 64, R = 146, and S 3 with N = 55, R = 120. We set the tensors corresponding to these fuzzy spaces explained in Section 3 (see Section 11 for the details of S 3 ) as the initial conditions at t = 0 of the equation of motion of the CTM shown in (42), and numerically obtained the solutions P abc (t) by the manner explained in Section 9. Then we performed the tensor-rank decompositions of P abc (t) for a number of representative values of t. Finally β(i, t) were determined by solving (36). Figure 9 plots the mean values, β(t) ≡ 1 R R i=1 β(i, t), in log-log plot. The gradients agree with −1 with the precisions down to the three decimal places, giving perfect agreement with (47).
As for B −1 µν (t), we have obtained the results shown in Figure 10. The left figure shows the time dependence of s max (4) defined in Section 8. This is expected to be proportional to B −1 µν , and therefore the data are fitted with (49). The agreement is rather nice with non-zero values of c R . This seems to contradict the supposed origin of c R , since S 1 does not have a curvature. The right figure shows the logarithmic derivative of the data obtained by subtracting the sequential data, and they are fitted with the corresponding derivative of (49). The results of the fitting are not nice but only barely acceptable.
The method above uses the short-time behavior of the virtual diffusion process, and is supposed to determine short distance structures of fuzzy spaces by measuring s max (i, j) between nearby points. If the method is fully reliable, one should be able to determine the time-dependence of the whole size of a homogeneous fuzzy space up to an allover factor by the local distance structure, because they should be proportional. On the other hand, we want to believe the validity of the continuum theory, because a fuzzy space with large N is made of many points and is expected to allow a continuum description. We seem to have a tension   (49)) with c R = 0.19, 0.24, respectively. Right: The logarithmic derivatives of the left data obtained by subtracting the sequential data. They are fitted with the logarithmic derivative of (49). This computation was not done for S 3 with N = 55, R = 120, since its size is too small to consider s max (4) as a reliable quantity (4 is well more than the half-size of S 3 ).
between the measuring method and the continuum theory.
To study the issue from a different angle, let us measure the whole size in a different manner using the lowest eigenvalue of a laplacian. The lowest non-zero eigenvalue of the minus of a laplacian is expected to be proportional to the inverse of the square size of a space. In our case, the inverse of the lowest non-zero eigenvalue of −K(i, j) in (37) is expected to behave in the same manner as B −1 µν . The left and the middle of Figure 11 plot the data of the inverse, and they are fitted with (49). The right one shows the logarithmic derivatives of the data obtained by subtracting the sequential data, and they are fitted with the logarithmic derivative of (49). They are in much better agreement with the continuum theory than Figure 10.
The two results above seem to conclude that the continuum theory is right, but our former method of measuring local distances for sizes is not fully reliable at least in our situation, while the latter method of using the lowest eigenvalue is. In fact, we have already pointed out the subtlety of the former method in Section 8. The former method would become more reliable and interesting for much larger fuzzy spaces with points existing more densely enough to validate the continuum description. There is also another possibility that the short-length dynamics is actually different from the global one. Note that the main contributors in the former method are the modes with large eigenvalues of −K in (37) because of short diffusion time, while the latter is the lowest one.
Another interesting thing in Figure 11 in comparison with Figure 10 is that the value of c R for S 1 has substantially decreased from the former method to the latter, while it keeps a more or less similar value for S 2 . This seems to suggest that a large portion of c R for S 1 comes from the small scale rather than the global scale. This would explain the presence of c R even for S 1 , which has no curvatures, in the following sense. As derived in [23], the right-hand side of the second equation in (45) actually contains the terms like β µµ β νν (∇ µ β)(∇ ν β)/β 2 and   β µν β µ ν (∇ µ ∇ ν β)/β, which can potentially contribute to c R . We have ignored these terms because of the homogeneity of the spaces, but the discreteness locally violates this assumption in short distances. Therefore, it does not seem obvious that these terms can really be ignored in short distances.
Though the agreement of the fitting in Figure 11 is really good especially in the small t region, there exist some deviations in the large t region, as can clearly be seen in the right figure. In fact, the continuum description cannot be expected to be right in this region, because the sizes of the fuzzy spaces are so large that the points exist sparsely on them, or, in other words, discreteness is macroscopic. To see this from a different viewpoint, we show the data from two fuzzy S 3 's with N = 55, R = 120 but with different damping factors. The damping factor used for the data S3N55R120 is e −l 2 /L 2 with L = 3, where l denotes the angular momenta of the modes. On the other hand, for S3N55R120LB, it is e −l(l+2)/L 2 with L = 4, where l(l + 2) comes from the eigenvalues of the Laplace-Beltrami operator on S 3 . It is clearly seen that the behaviors are qualitatively different from each other in the large t region. This implies that, while the small t region can well be described by the continuum theory irrespective of the damping factor, the large t region can be affected by small-distance details of fuzzy spaces. This seems to be consistent with the fact that the discrete structure is macroscopic in the large t region. It would be an interesting future problem how such deviations from the continuum can be described.

Generality of real symmetric three-way tensors
In section 3, we explained how to calculate the 3-way tensor P abc by using an example of 2sphere S 2 , and showed the realization of homogeneous fuzzy 2-sphere. In this section we show that this methodology can be applied to other fuzzy spaces by using some low dimensional manifolds with various topologies. The main point of this section is to show the generality of real symmetric three-way tensors by these demonstrations and consequently the generality of the CTM. One can explicitly see that real symmetric three-way tensors can in principle represent any spaces with free choices of dimensions and topologies. So the way such tensors realize spaces is essentially distinct from that in the other Euclidean tensor models [4,5,6,9], in which the dimensions of building simplicial blocks are supposed to be equivalent to the numbers of ways (the amount of indices) of tensors.
First let us summarize the procedure to construct a fuzzy space corresponding to a compact manifold M : 1. Take a coordinate x µ and a positive-definite metric g µν on the considering manifold M .
2. Prepare a set of real basis functions {f a (x)} on the manifold. It is convenient to impose an orthonormalization condition: for all combination a, b, where g = det(g µν ). Since the dimension of the function space is infinite in general, we have to choose a finite subset from the complete basis suited for a practical purpose. There are no general procedures for that, but in each individual case, there is a proper one. For example, let us suppose that basis functions f a (x) are taken to satisfy the Helmholtz equation, where ∆ is Laplace-Beltrami operator on the manifold M . Here m a plays the role as a "frequency" associated to each value a of the indices, and provides a natural way to choose a subset from the complete basis by {f a (x) | m 2 a ≤ Λ 2 } with some parameter Λ. This Λ determines the part of the basis which is considered, and effectively determines the value of N . The physical reason for considering such lower frequency modes than a cut-off is that we are interested in defining a space which is modified in a small scale but keeps its ordinary properties otherwise.
3. Define "regularized" functionsf a (x) from f a (x). There also exist a freedom in the way to regularize, but in the case f a (x) satisfies (51), a natural definition off a (x) is with some damping scale L. It is good to choose L Λ in general. Here, the damping factor can be another function damping with m 2 a or with a similar damping behavior. As discussed in Section 4, this regularization smoothens the cutoff and is important for locality of fuzzy spaces and good behavior of the virtual diffusion process.
4. Calculate P abc by using Then P abc defines a fuzzy space. It is expected that the tensor-rank decomposition of the tensor and connecting neighboring points will give a discretized counterpart of the manifold M .

Spheres
The square integrable functions on an n-dimensional sphere can be represented by a linear combination of n-dimensional (generalized) spherical harmonics. Therefore, let us take the generalized spherical harmonics as the set of the orthonormal basis functions on S n .

Circle S 1
We take for the coordinate on S 1 , θ ∈ [0, 2π] and √ g = 1. The set of basis functions is Figure 12: The homogeneous fuzzy 3-sphere with N = 55 and R = 120. The edges are drawn between the neighboring points by the criterion v i a v j a > 0.2, after the tensor-rank decomposition.
The notation like {sin nθ} is the abbreviation of {sin nθ| n ∈ N + }. Since these functions satisfy (51), we can use the procedure (52) to regularize the basis and the results arẽ with n ∈ N + . P abc can be calculated by and one can get homogeneous fuzzy circles, which look like polygons, from this P abc .
11.1.2 Three-dimensional sphere S 3 Figure 12 shows a homogeneous fuzzy 3-sphere obtained from the three-way tensor constructed from the above procedure. For this, we took n = 3, and l 3 was taken up to l 3 ≤ 4, which resulted in N = 55. The tensor-rank decomposition was carried out with R = 120, and the points have connections if v i a v j a > 0.2 in Figure 12. Though it is really hard to recognize this object as S 3 , the topological data analysis method discussed in Section 5, namely, persistent homology, is quite helpful. The analysis of Betti intervals is shown in Figure 13. The result tells the homology groups to be and dim(H 1,2,4 ) = 0. The homology groups agree with those of the 3-sphere, supporting our construction procedure of the fuzzy 3-sphere.

Line segments
One can also consider manifolds with boundaries. In the case of spaces with boundaries, one needs to set boundary conditions for its basis functions. There exist various choices such as Dirichlet, Neumann, their mixtures, and so on. In the following analysis, we simply consider the standard Dirichlet and Neumann boundary conditions.
Let us consider line segments. We take the coordinate of a line segment to be given by x ∈ [−π, π] and √ g = 1. In the case of Dirichlet boundary condition, one imposes f D a (±π) = 0 for all a, and then one finds two types of functions: sin(nx) and cos((n−1/2)x) where n ∈ N + . Thus the set of an orthonormal basis can be taken to be In the case of Neumann boundary condition, d dx f N a (x) x=±π = 0, there are three types of functions: the constant function, cos(nx) and sin((n − 1/2)x) with n ∈ N + . Then the set of an orthonormal basis can be taken to be The P abc can be computed from the regularized functionsf N,D a (x). The regularization factor can be taken for example to be exp(−k 2 /L 2 ) for a trigonometric function with frequency k. From these regularized functions, the three-way tensors can be obtained by where we are supposed to take a finite number of low-frequency modes. We have explicitly checked that fuzzy line segments are obtained from the tensor-rank decompositions of these tensors for both Dirichlet and Neumann boundary conditions.

Fiber bundles
In this subsection, let us construct more non-trivial fuzzy spaces, namely from fiber bundles.

Trivial bundles
Let us consider a manifold M which is isomorphic to the Cartesian product of two manifolds M 1 and M 2 . There exist the following relations: The basis functions f a (x) = f a 1 (x 1 )f a 2 (x 2 ) on M are normalized properly by the normalizations on M 1 and M 2 .

Möbius strip
The Möbius strip is an example of a nontrivial bundle, which is a bundle of a line segment over a circle. It can practically be built by considering a square and gluing a pair of opposite edges with a twist. From this we can easily find the conditions on the functions f (x, y) on a Möbius strip. We assume x, y ∈ [−π, π] and suppose that the edges at y = ±π are glued with a twist in the x direction. Then the condition on these edges gives the periodic boundary condition, The boundary condition on x = ±π can be freely chosen for instance from the Dirichlet boundary condition, or the Neumann boundary condition, for all y ∈ [−π, π]. Using the periodic condition (69) twice, one has with n, m ∈ N + . The regularized basis functionsf a (x, y) can be obtained by the procedure (52), and one obtains the three-way tensors by We have checked that connecting neighboring points by the result of the tensor-rank decompositions of the P produces discrete analogues of the Möbius strip for both Dirichlet and Neumann boundary conditions.

Klein bottle K 2
The Klein bottle is another nontrivial bundle of a circle over a circle and can be constructed by considering a square, gluing one pair of opposite edges, and gluing the other pair with a twist. This procedure tells us how to obtain the basis functions on a Klein bottle. We again take the coordinates (x, y) ∈ [−π, π] 2 and suppose that the glued edges with twisting correspond to those at y = ±π. Then the periodic boundary conditions are given by Using the condition (77) twice, one obtains So we see that f (x, y) can be expanded in a linear combination of exp(inx+imy/2) (n, m ∈ Z). Taking a real basis and requiring (77) with n, m ∈ N + . By the procedure explained before at (52), one can obtain regularized basis functionsf a (x) and the tensor P abc from them. Figure 14 shows the fuzzy Klein bottle with N = 49, the damping scale L = 3, and the tensor rank R = 49. Here N = 49 comes from setting the parameter Λ below (51) by Λ = 4. More explicitly, in the case Λ = 4, 49 is the summation of the numbers of the modes as N = 1 + 4 + 4 + 4 + 8 + 10 + 8 + 10 = 49, where the summands are ordered in the same way as in the expression (79). Note that the numbers of the combinations (n, m) ∈ N + × N + which satisfy n 2 + m 2 ≤ 4 2 and n 2 + (m − 1/2) 2 ≤ 4 2 are 8 and 10, respectively. The object in Figure 14 can be seen as a discretized two-dimensional closed surface with the structure of   self-intersection, which is the characteristics of Klein bottle. Figure 15 shows the Z 2 -coefficient Betti intervals for the fuzzy Klein bottle. This result shows that the homology groups of the fuzzy space are and we have also checked H n (K 2 , Z 2 ) = 0 at least for n = 3, 4, 5. Figure 16 shows the Z 3coefficient Betti intervals. This result also shows that the homology groups of the fuzzy space are and we have also checked H n (K 2 , Z 2 ) = 0 at least for n = 3, 4, 5. These homology groups agree with those of the ordinary continuous Klein bottle, supporting the validity of our construction of the fuzzy Klein bottle.
The contents of this section are limited to homogeneous fuzzy spaces, and we are successful at least in these cases. However, the construction procedure explained at the beginning of this section (also in Section 3) is not limited to homogeneous spaces, and it should be straightforward to construct inhomogeneous fuzzy spaces in a similar manner.

Summary and future prospects
The canonical tensor model (CTM) is a discrete model of gravity, which has a canonical conjugate pair of real symmetric three-way tensors as its dynamical variables. A question about the model was how to interpret the tensors as spacetimes. We have solved this question by using two well-known techniques in data analysis, namely the tensor-rank decomposition and persistent homology, and have formulated a mathematical procedure to extract topological and geometric properties from the real symmetric three-way tensors. We have also provided a systematic method to construct real symmetric three-way tensors corresponding to fuzzy spaces with any dimensions and topologies. We demonstrated these techniques by considering the real symmetric three-way tensors corresponding to homogeneous fuzzy S 1 , S 2 , and S 3 , solved the equations of motion of the CTM with these tensors as the initial conditions, and interpreted the time-dependent solutions as time-evolutions of geometric spaces. We have found that the results coincide with the expectation from the general relativistic system derived previously in a formal continuum limit of the CTM [23]. We have also explicitly constructed real symmetric three-way tensors for a variety of homogeneous fuzzy spaces with various dimensions and topologies, demonstrating the generality of the construction and extraction procedures, and hence of the CTM.
It is now apparent that the CTM is not limited to a particular dimension: The real symmetric three-way tensors can represent spaces of any dimension. This is a strong advantage of the CTM in relation to quantum gravity, because now generic spacetimes, including their dimensions, should emerge from the dynamics in the macroscopic limit rather than input parameters. This is in sharp contrast with the other Euclidian types of tensor models [4,5,6,9], in which the numbers of ways (indices) of tensors are directly related to the dimensions of simplicial building blocks of spaces.
Another important implication of this paper is that the formal continuum limit of the CTM discussed previously in [21,23] can actually be realized in large N limits. In these previous papers, the limit was formally put in by hand by performing an immediate replacement of discrete values of indices to continuum ones. In other words, the arguments were valid after the emergence of continuous macroscopic spacetimes, but did not tell anything about how they emerged. Though this paper is limited to the classical cases, we have explicitly shown that such limits can be realized by some large N cases by choosing appropriate initial conditions representing fuzzy spaces for the classical equation of motion. An obvious remaining problem here is how such initial conditions and classical trajectories are generated in the quantum framework, and we have a plausible hint for this: The physical wave function of the CTM has strong peaks at the tensor configurations invariant under Lie-groups with indefinite signatures [20,55]. As in the constructions of homogeneous spaces, such Lie-group symmetries can play vital roles in spacetime emergence.
Though this paper has introduced some interesting tools to interpret the dynamics of the tensors in the CTM as the dynamics of spacetimes, the applications are largely immature. This paper only dealt with the dynamics of zero modes in spaces, but for real physical interests, one has to deal with local dynamics in three dimensional spaces. Though there are no theoretical difficulties, there is a technical issue: Such studies require much larger fuzzy spaces, but the present performance of the tensor-rank decomposition is too slow. We are aware of high interest in studies in this direction of computer science, and hope that we are able to overcome this main technical difficulty in near future by incorporating recent developments. Other directly related future research directions would be trying to find and examine observables on the discrete geometry generated by the three-way tensors. Examples of such observables are already known from different discrete approaches, such as the spectral dimension [51], volume profile, Hausdorff dimension or the recently introduced quantum Ricci curvature [56,57]. It would also be interesting to see if this method can be extended to include matter fields to see the influence of matter on the dynamics of the fuzzy spaces. Furthermore, we have an interpretation for one of the tensors, but what would the geometric interpretation of the canonical conjugate be? It would also be interesting to see if the vectors from the tensor rank decomposition would have a quantum mechanical analogue.
The novel connection between gravity and data analysis shown in this paper stimulates some new kinds of, more speculative, questions. Can the Universe purely be described by data? How can one identify physically significant observables from random data? How do black holes appear in data? What is mass or energy in data? Is the equation of motion of the CTM useful in data analysis? We hope mutual communications of ideas in different fields stimulate new questions and studies to benefit them altogether. This is done in the latter part of the program. Each v i (i = 1, 2, . . . , R) is iteratively improved with the presence of the other vectors by minimizing the size of is kept fixed during the optimization of v i . This iteratively goes through i = 1, 2, . . . , R, forming one cycle. After every cycle, it is checked whether P abc v i a v i b v i c > 0 is satisfied by every i or not. If not, the v i which does not satisfy the condition is discarded. This cycle is repeated many times until the error ∆P abc cannot be reduced or becomes smaller than a criterion.
In our application, it is observed that the above condition P abc v i a v i b v i c > 0 tends to avoid rough tensor-rank decompositions containing mutual cancellations among large v i 's 21 , and gives a decomposition with v i of nearly equal sizes. This is useful in our application, because we are considering homogeneous fuzzy spaces, in which all the points should be more or less uniformly weighted.
Finally let us explain the actual minimization method used as a subroutine in the above procedure. The subroutine minimizes the size of for a givenP by optimizing v. By taking the square, (P abc − v a v b v c )(P abc − v a v b v c ), this is to find a minimum of Though this is a well-defined problem, it is not easy to obtain a global minimum, and we restrict ourselves to finding a local one. This limitation is a disadvantage which cannot be underestimated in general, but, since this subroutine is called many times, it is also unlikely that v stays in a bad local minimum through the whole process. The condition for a local minimum of (91) is given by the vanishing of its first derivative, and the absence of negative eigenvalues of its second derivative matrix (Hessian) given by Rather than trying to directly solve the above problem which has strong non-linearity, let us consider a simpler form than (91), 21 See for example [37,42] for more details about this numerically (and theoretically) serious problem.
The minimization of this is also well-defined. The reason for considering (94) rather than (91) becomes evident in due course. In the same way as above, the conditions for a local minimum of (94) are given by and the non-negativity (in the same meaning as above) of 2w a w b + (w 2 )δ ab − 2P abc w c .
Comparing the two problems, one can see that a local minimum of the latter gives one of the former by doing a rescaling v 2 v a = w a . Here it is important that the non-negativity of (96) readily implies that of (93), because the difference 2v 2 v a v b is non-negative.
A way to obtain a local minimum of (94) is to apply the steepest descent method. By taking the first derivative of (94) and choosing a step size of γ/w 2 with γ > 0, one obtains a sequence, A convergent vector of the sequence gives a local minimum of (94).
An advantage of considering (97) from (94) rather than what can be obtained from (91) is that (97) is more controllable. The second term of (97) is bounded for any w, and one can easily prove that, for 0 < γ < 2, the sequence does not diverge. Therefore, the worst behavior is a bounded non-convergent sequence, which would be changed to a convergent one by an appropriate choice of γ, typically by making it smaller. In our application, however, the simplest choice γ = 1 suffices 22 , where only the second term needs to be computed on the righthand side of (97).