Causal Capsules and Tensor Autoencoders
Year: 2021

Abstract:
This paper introduces a set of neural network architectures for forward and inverse causal inference that are consistent with capsule theory and implement multilinear (tensor) factor analysis methods. Forward causal inference is addressed with a causal autoencoder-decoder architecture composed of a set of causal capsules that estimate the latent variables representing the constituent factor of data formation, and a tensor-autoencoder that governs the latent variable interaction. A recurrent non-linear causal capsule chain that employs kernel activations computes the optimal linearized subspace for every causal factor, and implements the kernel multilinear principal component analysis or the kernel multilinear independent component analysis. For distributed computation, we break the chain links and each causal representation is computed separately, shuttling causal information between capsules. The causal factor representations may be computed efficiently by restructuring the input into a hierarchy of parts with a set of part-based causal capsules that are “glommed” together to create a part-based hierarchy of causal capsules. Inverse causal inference, the estimation of causes of effects, is addressed with a multilinear projection architecture that inverts the estimated forward causal model and employs a set of observations to constrain the solution set rendering the problem well-posed.

heading: 1 INTRODUCTION
text: Neural networks are being employed increasingly in high-stakes application areas, such as face recognition [Taigman et al. (2014); Huang (2012); Sun et al. (2013); Chen et al. (2015); Xiong et al. (2016)], and medical technologies [Kermany et al. (2018); Madani et al. (2018); Topol (2019)]. Developing a set of neural network architectures that are causally explainable is important in developing a trustworthy machine learning, where "A causes B" means "the effect of A is B", a measurable and experimentally repeatable quantity [Holland (1986)].
Forward causal inference models the mechanism of data formation, and estimates the effects of interventions [Pearl (2000); Imbens & Rubin (2015); Spirtes et al. (2000); Vasilescu et al. (2021); Vasilescu & Terzopoulos (2002a; 2005; 2004)]. Unlike, conventional statistics and conventional machine learning that model the observed data distribution, and make predictions about a variable that has been co-observed with another. Inverse causal inference estimates the causes of effects given an estimated forward causal model that is inverted subject to a set of observations that constrain the solution set [Vasilescu (2011); Vasilescu & Terzopoulos (2007)].
There are two conceptual frameworks for causal inference: DAGs or path analysis and potentialoutcome. Donald Rubin and his collaborators have advocated the potential outcome approach which framed causal inference as a missing data problem [Imbens (2020)]. Judea Pearl has been advocating do-calculus – a directed acyclic graph approach as a mathematical language that he has unified with structural equations and counterfactuals [Bollen & Pearl (2013)]. Judea Pearl’s causation ladder [Pearl (2000)] provides a way of thinking about causal discovery, causal reasoning, and decision making. Pearl & Bareinboim (2014); Bareinboim & Pearl (2016) have parameterized the differences between experimental and observational studies based on possible sources of error.
Tensor data analysis is a type of structural equation modeling that has been employed to perform dimensionality reduction, to develop regression models, and to model cause-and-effect, Fig. 1. Tensor factor analysis has been employed in psychometrics [Tucker (1966); Harshman (1970); Carroll & Chang (1970); Bentler & Lee (1979); Kroonenberg & de Leeuw (1980)], econometrics [Kapteyn et al. (1986); Magnus & Neudecker (1988)], chemometrics [Bro (1997); Acar et al. (2014)], signal processing [de Lathauwer (1997; 2008); Cichocki et al. (2009)], computer vision [Vasilescu & Terzopoulos (2002b); Wang & Ahuja (2003)], computer graphics [Vasilescu (2002); Davis & Gao (2003); Vasilescu & Terzopoulos (2004); Vlasic et al. (2005); Hsu et al. (2005)], and machine learning [Vasilescu (2009); Vasilescu & Terzopoulos (2005)]. In machine learning, tensor methods have been effectively employed to reparameterize neural networks. Neural network weights have been organized
into “data tensors” and dimensionally reduced in order to achieve greater computational efficiency [Lebedev et al. (2014); Novikov et al. (2015); Kim et al. (2015); Khrulkov (2020); Onu et al. (2020); Iwen et al. (2021)]. Tensor methods have also been applied in regression analysis [Kolda et al. (2005); Chu & Ghahramani (2009); Tang et al. (2013); Anandkumar et al. (2014); Kossaifi et al. (2017); Wang et al. (2017); Benesty et al. (2021); Vendrow et al. (2021)].
This paper introduces a set of causal capsule architectures for forward and inverse causal inference that implement tensor factor analysis operations. These architectures are consistent with capsule theory proposed by Geoffrey Hinton and his collaborators [Hinton et al. (2011); Sabour et al. (2017); Hinton (2021)].
Forward causal inference, the estimation of effects of causes, is performed with a causal autoencoder architecture that consists of several causal capsules that compute the causal factor representations, and a tensor-autoencoder that governs the causal factor interaction, Fig 2. A causal capsule is formed from a set of constrained “cluster”-based autoencoders1that transform the basis vectors spanning the “cluster” subspace, such that a causal factor representation is invariant of the “cluster” membership, i.e., invariant to all the other causal factors [Vasilescu (2009)]. A tensor autoencoder is an autoencoder with a vectorized tensor code formed from the multilinear (tensor) product of factor representations.
A recurrent non-linear causal capsule chain that employs kernel activations computes the optimal linearized subspace for every causal factor, and implements the kernel multilinear principal component analysis or the kernel multilinear independent component analysis, Fig 3. For distributed computation, we break the chain links andeach causal representation is computed separately, shuttling causal information between capsules.
For a scalable architecture, causal representations for an object whole can be computed efficiently by parts, Fig 4. As Hinton (2021) has also indicated, a part-based causal capsule architectures may also be “glommed” together to analyze a hierarchy of data columns [Vasilescu et al. (2021); Vasilescu & Kim (2019); de Lathauwer (2008)]. The hierarchical neural network architecture is a compositional
1In the context of multifactor data analysis, a cluster is a set of observations for which all factors are fixed except one. Data belonging to the same cluster may not form a cluster in Euclidean space and not easily identifiable by an EM algorithm [ Dempster et al. (1977)].
hierarchical computation of causal factor representation and implements the Incremental M -mode Block SVD algorithm. 2
Algorithm 1 M -mode SVD algorithm. Input the data tensor D ∈ CI0×···×IM .
1. For m := 0, . . . ,M , Let Um be the left orthonormal matrix of [UmSmVTm] := svd(D[m]) a 2. Set Z := D ×0 U0T ×1 U1T · · · ×m Um
T · · · ×M UM T. Output mode matrices U0,U1 . . . ,UM and the core tensor Z .
aThe computation of Um in the SVD D[m] = UmΣVm
T can be performed efficiently, depending on which dimension of D[m] is smaller, by decomposing either D[m]D[m]T = UmΣ2UmT (note that VmT = Σ+UmTD[m]) or by decomposing D[m] TD[m] = VmΣ 2Vm
T and then computing Um = D[m]VmΣ +.
Inverse causal inference is performed with a multilinear projection architecture [Vasilescu & Terzopoulos (2007); Vasilescu (2009)] that is performed by inverting an estimated forward model subject to data constraints. Fig. 5.
The architectures are derived based on two mathematical principles: (i)linear autoencodersdecoders weights are the principal component analysis basis vectors, sec. 2, (ii) the object-whole representation can be derived bottom-up in closed form from a part-based hierarchical causal factor representation, sec 3.2.
After reviewing the mathematical foundations of our work in the next section, we discuss forward causal models and depict their neural network arcitectures in Section 3 and discuss inverse causal inference and depict the multilinear projection neural network architectures in Section 4. Section 5 concludes the paper.

heading: 2 LINEAR AUTOENCODER AND LINEAR PCA
text: An autoencoder-decoder that minimizes the reconstruction loss function,
l = I∑ i=1 ‖di −Bci‖+ λ‖BTB− I‖ (1)
and has a linear decoder learns a set of weights, br that are identical to the PCA basis vectors when the weights of each neuron, cr are computed sequentially. An autoencoder is implemented with a cascade of Hebb neurons [Hebb (1949)]. The contribution of each neuron, c1, . . . , cr, are the PCA sequentially computed and subtracted from a centered training data set, and the difference is driven through the next Hebb neuron, cr+1 [Sejnowski et al. (1989); Sanger (1989); Rumelhart et al. (1986); Ackley et al. (1985); Oja (1982)]. The weights of a Hebb neuron, cr, are updated by
∆br(t+ 1) = η
( d−
r∑ ir=1 bir(t)cir(t)
) cr(t) = η ( d−
r∑ ir=1 bir(t)b T ir(t)d
) dTbr(t),
br(t+ 1) = (br(t) + ∆br(t+ 1))
‖br(t) + ∆br(t+ 1)‖ where d ∈ CI0 is a vectorized centered observation with I0 measurements, η is the learning rate, br are the autoencoder weights of the r neuron, cr is the activation, and t is the time iteration. Backpropagation [LeCun et al. (1988; 2012)] is equivalent to performing PCA gradient descent [Jolliffe (1986)].

heading: 3 CAUSAL INFERENCE
text: Throughout this article, we will denote scalars by lower case italic letters (a, b, ...), vectors by bold lower case letters (a,b, ...), matrices by bold uppercase letters (A,B, ...), and higher-order tensors by bold uppercase calligraphic letters (A,B, ...). Index upper bounds are denoted by italic uppercase letters (i.e., 1 ≤ a ≤ A or 1 ≤ i ≤ I). The zero matrix is denoted by 0, and the identity matrix is denoted by I. The TensorFaces paper [Vasilescu & Terzopoulos (2002a) is a gentle introduction to tensor factor analysis, Kolda and Bader [Kolda & Bader (2009) is a nice survey of tensor methods and references [Vasilescu (2009); de Lathauwer (1997); Bro (1997) provide an in depth treatment of tensor factor analysis.
2By comparison, a hierarchical Tucker is a resource efficient hierarchical computational scheme that employs a hierarchical re-balancing of the modes trick in which one flattens a data tensor in multiple modes at the same time to avoid computing SVDs of skinny matrices [Hackbusch & Kühn (2009); Grasedyck (2010); Perros et al. (2015)].
Algorithm 2 Kernel Multilinear PCA/ICA (K-MPCA/MICA) algorithm. Input the data tensor D ∈ CI0×···×IM , where mode m = 0 is the measurement mode, and the desired ranks R̃1, . . . , R̃M .
1. For m := 1, . . . ,M , Compute the elements of the mode-m covariance matrix, for j, k := 1, . . . , Im, as follows: [D[m]D[m] T]jk :=
I1∑ i1=1 ... Im−1∑ im-1=1 Im+1∑ im+1=1 ... IM∑ iM=1 K(di1...im-1 j im+1...iM ,di1...im-1 k im+1...iM).
 For K-MPCA: Set Um to the left matrix of the SVD of D[m]D[m]T = UmΣ2UmT Truncate to R̃m columns Um ∈ CIm×R̃m .
For K-MICA: Compute Um := Cm ∈ CIm×R̃m based on [Vasilescu & Terzopoulos (2005)]. The initial SVD truncates to R̃m.
2. Set T := D ×1 U+1 · · · ×m U+m · · · ×M U + M . 3. Local optimization via alternating least squares: Iterate for n := 1, . . . , N
For m := 1, . . . ,M , Set Xm := D ×1 U+1 · · · ×m−1 U + m−1 ×m+1 U + m+1 · · · ×M U + M .
Set Um to the R̃m leading left-singular vectors of the SVD of Xm,[m]a. Set T := XM ×M U+M .
until convergence. Output the converged extended core tensor T ∈ CI0×R̃1×···×R̃M and causal factor mode matrices U1, . . . ,UM .

heading: 3.1 FORWARD CAUSAL INFERENCE
text: Forward causal inference frames questions in terms of interventions. What if the causal factor c were changed by one unit, how much would the observed measurements, d, be expected to change?
For modeling individual level-effects rather than the average effects of causes, observations are acquired by systematically varying each causal factor while holding the rest of the causal factors fixed.
Within the tensor mathematical framework, a M -way array or “data-tensor”, D ∈ CI0×I1···×Im···×IM contains a collection of vectorized and centered observations,3 di1...im...iM ∈ RI0 that are the result of M causal factors. The m causal factor (1 ≤ m ≤ M ) takes one of Im values that are indexed by im, 1 ≤ im ≤ Im. An observation that is result of the confluence M causal factors is modeled by a multilinear tensor equation with multimode latent variables,
di1,...,iM = T ×1 (r T i1 + Ti1) · · · ×M (r T iM + TiM) + i1,...,iM
where T = Z×0U0 is the extended core which modulates the interaction between the latent variables, ri1 . . . rim . . . riM , that represent the causal factors, and im ∈ N (0,Σm) is an IID Gaussian noise.
3Reference [(Vasilescu, 2009, Appendix A)] evaluates some of the arguments found in highly cited publications in favor of treating an image as a matrix (tensor) rather than a vector. While technically speaking, it is not incorrect to treat an image as a matrix in linear/tensor algebra, most arguments do not stand up to analytical scrutiny, and it is preferable to vectorize an image and treat it as a single observation rather than a collection of independent column/row observations.
The M -mode SVD, Alg. 1 and Alg. 2, and their neural network architecture counterparts, Fig. 2, Fig. 3 may be employed to represent data in terms of their causal factors by minimizing the reconstruction loss function,
l = ‖D − T ×1 U1 · · · ×m Um · · · ×M UM‖+ M∑
m=1
λm‖UmUTm − I‖, (2)
where T is the extended core and the mode matrices, Um, spans the m causal factor representation. Each mode matrix, Um is computed using alternating least squares where a set of M least squares are computed by moving the mode matrices U1, . . . ,Um-1,Um+1, . . . ,UM to the knowns side of the equation, setting Xm := D ×1 · · · ×m-1 UTm-1 ×m+1 UTm+1 · · · ×M UTM for distributed computation or setting Xm := (Xm-1 ×m-1 UTm-1)×m Um) for sequential computation, and optimizing the loss function
l = ‖Xm − T ×m Um‖+ λ‖UmUTm − I‖ = ‖Xm[m] −UmT[m]‖+ λ‖UmUTm − I‖, (3) where Xm := D ×1 · · · ×m-1 UTm-1 ×m+1 UTm+1 · · · ×M UTM - distributed computation (4)
= (Xm-1 ×m-1 UTm-1)×m Um = T ×m Um - sequential computation (5)
The mode matrix Um is set to the subspace of the matrixized Xm, UmSmVTm := svd(Xm[m])4. Figure 3(e) displays a recurrent causal chain that unrolls the for-loop from step 3, Alg 2 and sequentially computes the mode matrices by performing an SVD on the a matrixized Xm computed from eq.( 5). For a distributed computation, M different Xm are computed according to eq.( 4), where mode matrices are shuttled between the different threads.

heading: 3.2 DERIVATION: INVARIANCE AND HIERARCHY OF CAUSAL CAPSULES
text: In this section, takes advantages of the principle that an SVD can be computed from its parts. On that basis, we derive a causal factor representation that is statistical invariant to “cluster” membership (i.e., all other causal factors). On that basis, we provide a scalable architecture by deriving a compositional bottom-up computation of an object whole representation. Thus, the naive and direct implementation of the M-mode SVD is replaced with a compositional hierachical part-based distributed architecture.
Computing the mode matrices, Um, may be viewed as equivalent to computing a set of mutually constrained, cluster-based PCAs. When dealing with data that can be separated into clusters, the standard machine learning approach is to compute a separate PCA. When data from different clusters are generated by the same underlying process (e.g., facial images of the same people under different viewing conditions), the underlying data can be concatenated in the measurement mode and the common causal factor can be modeled by one PCA.5
Thus, we define a constrained, cluster-based PCA as the computation of a set of PCA basis vectors that are rotated such that the latent representation is constrained to be the invariant of the cluster.
MPCA performs M constrained, cluster-based PCAs, since the computation of the mode-m matrix Um, which involves a mode-m data tensor flattening and subsequent SVD, is equivalent to performing a constrained, cluster-based PCA; i.e., data cluster concatenation followed by an SVD. This is self evident when employing our modified datum-centric flattening operator, Fig. 7.
In the context of our multifactor data analysis, we define a cluster as a set of observations for which all factors are fixed except one, the m factor. Note that there are Nm = I1I2 . . . Im−1Im+1 . . . IM possible clusters and the data in each cluster varies with the same causal mode.6 Thus, the data across different clusters share one of the underlying causal factors. The constrained, cluster-based PCA concatenates the clusters in the measurement mode and analyzes the data with a linear model, such as PCA or ICA [Bartlett et al. (2002); Common (1994); Lathauwer et al. (1995); De Lathauwer et al. (1996); Anandkumar et al. (2014).
To see this, let Di1...im−1im+1...iM ∈ CI0×1×1···×1×Im×1···×1 denote a subtensor of D that is obtained by fixing all modes except causal factor mode m and mode (the measurement mode). Matrixizing this subtensor in the measurement mode 0, we obtain Di1...im−1im+1...iM [0] ∈ C
I0×Im . This data matrix comprises a cluster of data obtained by varying the m causal factor, to which one can traditionally
4See Alg. 1, footnote a 5The active appearance model concatenated two different measurements, facial feature locations and texture, to compute a person representation for a particular viewpoint invariant of measurement. More generally, one can concatenate the measurements of a person from different viewpoint clusters to compute a person representation that is invariant of the measurement and the viewpoint [Cootes et al. (2001); Si et al. (2013)].
6Observations in the cluster may not be in Euclidean proximity in the measurement space. Consequently, a cluster may not be easily identified through a standard EM algorithm.
apply PCA. Since there are Nm = I1I2 . . . Im−1Im+1 . . . IM possible clusters that share the same underlying space associated with the cmth factor, the data can be concatenated and PCA performed in order to extract the same representation for the m factor regardless of the cluster. Now, consider the MPCA computation of mode matrix Um, Fig. 3(a), which can be written in terms of matrixized subtensors as
Dm T = 
D1...11...1[m] T ... DI1...11...1[m] T
... DI1...Im−1Im+1...IM [m] T
 = VmΣmUm T. (6)
Clearly, this is equivalent to computing a set of Nm = I1I2 . . . Im−1Im+1 . . . IM cluster-based PCAs concurrently by combining them into a single statistical model and representing the underlying causal factor m common to the clusters. Thus, rather than computing a separate linear PCA model for each cluster, MPCA concatenates the clusters into a single statistical model and computes a representation (coefficient vector) for mode m that is invariant relative to the other causal factor modes 1, ..., (m− 1), (m+ 1), ...,M . Thus, MPCA is a multilinear, constrained, cluster-based PCA.
To clarify the relationship, let us number each of the matrices Di1...im−1im+1...iM [m] = D (n) m with a parenthetical superscript 1 ≤ n = 1 + ∑M k=1,k 6=m(in − 1) ∏k−1 l=1,l 6=m Il ≤ Nm.
Let each of the linear SVDs be
D(n)m = U (n) m Σ (n) m U (n) 0 T (7) (8)
D[m] = [ U (1) m Σ (1) m . . . U (NM ) m Σ (Nm) m ] ︸ ︷︷ ︸
SVD
diag([ U (1) 0 . . . U (Nm) 0 ])T, (9)
= UmΣmW T m diag([ U (1) 0 . . . U (Nm) 0 ])T, (10)
= UmΣm[ U (1) 0 W (1) m . . . U (Nm) m W (Nm) m ]T (11)
where diag(·) denotes a diagonal matrix whose elements are each of the elements of its vector argument. The mode matrix U(nm)0 is the measurement matrix that contains the eigenvectors that span the observed data in cluster nm, 1 ≤ nm ≤ Nm. MPCA can be thought as computing a rotation matrix, Wm, that contains a set of blocks W(n)m along the diagonal that transform the PCA cluster eigenvectors, U(nm)0 , such that the mode matrix Um is the same regardless of cluster membership, eqs.(9-11). The constrained “cluster”-based PCAs may also be implemented with a set of concurrent “cluster”-based PCAs.
Object wholes appear to have have a hierarchy of perceptual parts. Part-based causal capsule architectures may be “glommed” together to create a hierarchy of part-based causal capsules that analyze a hierarchy of data columns [Vasilescu et al. (2021); Vasilescu & Kim (2019); de Lathauwer (2008)]. This hierarchical architecture implements the incremental hierarchical multilinear (tensor) block decomposition algorithm.
Causal factors of object wholes may be computed efficiently from their parts, Fig 4. The matrixized data tensor may be organized into part “clusters” by applying a permutation
D ×T m P ⇔ PD[m]T (12) where P is a permutation matrix. The resulting hierarchical architecture implements the Incremental M-mode Block SVD algorithm. The Incremental M-mode Block SVD is a generalized hierarchical part-based decomposition that computes an exact global decomposition and is suitable for streaming data [Vasilescu et al. (2021).

heading: 3.3 NONLINEAR CAUSAL CAPSULES AND KERNEL MPCA/KERNEL MICA
text: An autoencoder with a non-linear activation function represents an observation with
di = fd(Bd fe(Bedi + ae)︸ ︷︷ ︸ ci +ad), (13)
where fe, fd are the encoder, decoder activation functions, and Be,Bd are the encoder, decoder weights respectively. Kernel PCA (KPCA) is often given as an example of a “true” nonlinear model. KPCA first applies a nonlinear transformation to the data and then it performs a linear decomposition. Thus, KPCA derives its nonlinearity from its preprocessing step. Other nonlinear methods include nonlinear PCA (NLPCA) [Kramer (1991), as well as kernel PCA (KPCA) [Schölkoph et al. (1998) and kernel LDA (KLDA) [Yang (2002) methods in which kernel functions that satisfy Mercer’s theorem correspond to inner products in infinite-dimensional space. An alternative approach is to apply linear models to nonlinear problems through the “kernel trick”, specifically the kernel PCA [Schölkoph et al. (1998) and kernel ICA [Yang et al. (2005) techniques.7 Kernel PCA/ICA are
7The so-called “kernel trick” maps the original non-linear measurements into a higher-dimensional space, where a linear classifier is subsequently used; this makes a linear classification in the new space equivalent to non-linear classification in the original space. This is done using Mercer’s theorem, which states that any continuous, symmetric, positive semi-definite kernel K(u,v) can be expressed as an inner product in a highdimensional space. Wherever an inner product between two vectors is used, it is replaced with a kernel of the
a
b
cFigure 4: (a) Learning levels of abstractions bottom-up with a hierarchy causal part-based capsules. Each causal capsule in Fig 3 (f) can be replaced with a hierarchy of parts and wholes. The constrained cluster-based PCA (b), may be computed with a hierarchy of autoencoders, as derived in eq. (6)-(10), and depicted in Fig 3.
nonlinear versions of their conventional linear counterparts, but of course they are not multimodal factor models.
The kernel trick can also be applied to our multilinear, multifactor PCA/ICA models to further nonlinearize them, thus enabling them to deal with arbitrarily nonlinear data, Alg. 2.
To accomplish this, recall that the computation of the mode-m covariance matrix D[m]D[m]T involves inner products dTi1...im-1 j im+1...iMdi2...im-1 k im+1...iM between pairs of images in the image data tensor D associated with causal factor mode m, for m = 1, . . . ,M . We replace the inner products with a generalized distance measure between images, K(di1...im−1 j im+1...iM ,di2...im−1 k im+1...iM ), where K(·, ·) is a suitable kernel function (Table 1), which corresponds to an inner product in some expanded feature space. This generalization naturally leads us to a Kernel Multilinear PCA (K-MPCA) Algorithm, where the covariance step computation in in Algorithm 1 is replaced by
[D[m]D[m] T] jk := I1∑ i1=1 · · · Im−1∑ im−1=1 Im+1∑ im+1=1 · · · IM∑ iM=1 K(di1...im−1 j im+1...iM ,di1...im−1 k im+1...iM ).
Similarly, a Kernel Multilinear ICA (K-MICA) Algorithm results from making the same generalization in the MICA algorithm [Vasilescu & Terzopoulos (2005). Algorithm 2 specifies both K-MPCA and K-MICA. Figure 3(d) unrolls the for-loop in step 3 of Alg. 2.

heading: 4 INVERSE CAUSAL INFERENCE:MULTILINEAR AND MULTIPLE LINEAR PROJECTIONS
text: Inverse causal inference estimates the causes of effects, and addresses the why question. Inverse problems often violate one of the conditions of a well-posed problem, and Donald Rubin has referred to the "why" question as “cocktail party chatter” [Gelman & Imbens (2013). For a problem to be well-posed a solution must exit, it must be unique and the solution’s behaviour ought to change continuously with the initial conditions [Hadamard (1952). Often, there are multiple combinations of same causal factors that have the same potential outcome. In imaging, these types of outcomes are known as visual illusions.
Therefore, inverse causal inference is the estimation of causes of effects given an estimated forward causal model that is inverted subject to a set of observations that constrain the solution set and render the problem well-posed [Vasilescu et al. (2021)].Similar to reverse causal inference [Gelman & Imbens (2013)], inverse causal inference may be employed as a model checking mechanism and motivation for forward inference question.
vectors. Thus, a linear algorithm is easily transformed into a nonlinear algorithm. This trick has been applied to numerous algorithms in machine learning and statistics.
Multilinear projection simultaneously projects one or more unlabeled test images that are not part of the training data set into multiple constituent causal factor spaces associated with data formation, in order to infer the mode labels:
CP or M -mode SVD(T +x×Txdtest) ≈ r1...◦rm...◦rM◦rE. Topologically the multilinear projection architecture, Fig. 5, is an inverted M-mode SVD architecture. When the dimensionality of vec(R) is larger than the number of measurements in d, then the system of equations is under determined. There are three possible solutions – dimensionality reduction of the mode matrices, and modeling the mechanism of data formation by multiple linear or tensor models. Instead of performing a multilinear projection, [Vasilescu & Terzopoulos (2002b) perform a set of linear projections.

heading: 5 CONCLUSION
text: This paper introduces deep causal learning architectures that implement tensor factor analysis operations and model the mechanism of data formation. The tensor factor analysis methods, the M -mode SVD, the Kernel MPCA/MICA, and the associated causal capsules architectures transform the “cluster” eigenvectors such that the constituent causal representations are invariant of the cluster membership, i.e., invariant of other causal factors of data formation. Causal representation may be computed efficiently by “glomming” together a hierarchy of part-based causal capsules. The hierarchical part-based causal architecture implements the compositional hierarchical tensor factorization, the Incremental M-mode Block SVD. Each part-based capsule analyzes a data column from a hierarchy of data columns [Vasilescu et al. (2021). Inverse causal inference, the estimation of causes of effects, is accomplished with a multilinear projection algorithm. The neural architecture that implements the multilinear projection is an inverted M-mode SVD architecture. Tensor causal factor analysis and their associates neural networks have properties consistent with the capsule theory. Tensor causal factor analysis has been applied on real and synthetic data in many domains, including face recognition where the approach is known as TensorFaces and computer graphics where a set of TensorTextures are synthesized for arbitrary geometries.

heading: A MATHEMATICAL BACKGROUND
text: Throughout this article, we will denote scalars by lower case italic letters (a, b, ...), vectors by bold lower case letters (a,b, ...), matrices by bold uppercase letters (A,B, ...), and higher-order tensors by bold uppercase calligraphic letters (A,B, ...). Index upper bounds are denoted by italic uppercase letters (i.e., 1 ≤ a ≤ A or 1 ≤ i ≤ I). The zero matrix is denoted by 0, and the identity matrix is denoted by I. The TensorFaces paper [Vasilescu & Terzopoulos (2002a) is a gentle introduction to tensor factor analysis, [Kolda & Bader (2009) is a great survey of tensor methods and references [Vasilescu (2009); de Lathauwer (1997); Bro (1997) provide an in depth treatment of tensor factor analysis.
A.1 PCA COMPUTATION WITH LINEAR AUTOENCODER An autoencoder-decoder that minimizes the reconstruction loss function for a set of observations, di ∈ CI0 ,
l = I∑ i=1 ‖di −Bci‖+ λ‖BTB− I‖, (14)
and has a linear decoder learns a set of weights, bi0,r, that are identical to the elements of the PCA basis matrix, B ∈ CI0×R, when the weights of each neuron are computed sequentially, Fig. 6. An autoencoder is implemented with a cascade of Hebb neurons [Hebb (1949). The contribution of each neuron, c1, . . . , cr, is sequentially computed and subtracted from a centered training data set, and the difference is driven through the next Hebb neuron, cr+1 [Sejnowski et al. (1989); Sanger (1989); Rumelhart et al. (1986); Ackley et al. (1985); Oja (1982).
The weights of a Hebb neuron, cr, are updated by
∆br(t+ 1) = η
( d−
r∑ ir=1 bir(t)cir(t)
) cr(t) (15)
= η ( d−
r∑ ir=1 bir(t)b T ir(t)d
) dTbr(t),
br(t+ 1) = (br(t) + ∆br(t+ 1))
‖br(t) + ∆br(t+ 1)‖
where d ∈ CI0 is a vectorized centered observation with I0 measurements, η is the learning rate, br are the autoencoder weights of the r neuron, cr is the activation, and t is the time iteration. Back-propagation [LeCun et al. (1988; 2012) performs PCA gradient descent [Jolliffe (1986).
A.2 RELEVANT TENSOR ALGEBRA Briefly, the natural generalization of matrices (i.e., linear operators defined over a vector space), tensors define multilinear operators over a set of vector spaces. A “data tensor” denotes an M -way data array.
Definition 1 (Tensor) Tensors are multilinear mappings over a set of vector spaces, CIm , 1 ≤ m ≤ M , to a range vector space CI0 :
A : { CI1 × CI2 × · · · × CIM } 7→ CI0 . (16)
The order of tensor A ∈ CI0×I1×···×IM is M + 1. An element of A is denoted as Ai0i1...im...iM or ai0i1...im...iM , where 1 ≤ im ≤ Im.
The mode-m vectors of an M -order tensor A ∈ CI0×I1×···×IM are the Im-dimensional vectors obtained from A by varying index im while keeping the other indices fixed. In tensor terminology, column vectors are the mode-0 vectors and row vectors as mode-1 vectors. The mode-m vectors of a tensor are also known as fibers. The mode-m vectors are the column vectors of matrix A[m] that results from matrixizing (a.k.a. flattening) the tensor A.
Definition 2 (Mode-m Matrixizing) The mode-m matrixizing of tensor A ∈ CI0×I1×...IM is defined as the matrix A[m] ∈ CIm×(I0...Im−1Im+1...IM ). As the parenthetical ordering indicates, the
Algorithm 3 M -mode SVD algorithm. Input the data tensor D ∈ CI0×···×IM .
1. For m := 0, . . . ,M , Let Um be the left orthonormal matrix of [UmSmVTm] := svd(D[m]) a
2. Set Z := D ×0 U0T ×1 U1T · · · ×m UmT...×M UM T. Output mode matrices U0,U1 . . . ,UM , and the core tensor Z .
aThe computation of Um in the SVD D[m] = UmΣVmT can be performed efficiently, depending on which dimension of D[m] is smaller, by decomposing either D[m]D[m]T = UmΣ2UmT (note that VmT = Σ+Um TD[m]) or by decomposing D[m]TD[m] = VmΣ2VmT and then computing Um = D[m]VmΣ+.
mode-m column vectors are arranged by sweeping all the other mode indices through their ranges, with smaller mode indexes varying more rapidly than larger ones; thus,
[A[m]]jk= ai1...im...iM , where (17)
j = im and k = 1 + M∑ n=0 n 6=m (in − 1) n−1∏ l=0 l 6=m Il.
A generalization of the product of two matrices is the product of a tensor and a matrix [de Lathauwer (1997); Carroll et al. (1980).
Definition 3 (Mode-m Product, ×m) The mode-m product of a tensor A ∈ CI1×I2×···×Im×···×IM and a matrix B ∈ CJm×Im , denoted by A ×m B, is a tensor of dimensionality CI1×···×Im−1×Jm×Im+1×···×IM whose entries are computed by
[A×mB]i1...im−1jmim+1...iM= ∑ im ai1...im−1imim+1...iM bjmim ,
C = A×m B. matrixize
tensorize C[m] = BA[m].
The M -mode SVD, Alg. 1 [Vasilescu & Terzopoulos (2002a) is a “generalization” of the conventional matrix (i.e., 2-mode) SVD which may be written in tensor notation as
D = U0SU1 T ⇔ D = S×0 U0 ×1 U1
The M -mode SVD orthogonalizes the M spaces and decomposes a tensor as the mode-m product, denoted ×m , of M -orthonormal mode matrices, and a core tensor Z
D=Z ×0 U0 · · · ×m Um · · · ×M UM . (18) D[m] = UmZ[m] (UM · · · ⊗Um+1 ⊗m-1 U · · · ⊗U0) T, (19) vec(D)=(UM · · · ⊗Um+1 ⊗Um-1 · · · ⊗U0) vec(Z). (20)
The latter two equations express the decomposition in matrix form and in terms of vec operators.
A.3 COMPOSITIONAL HIERARCHICAL BLOCK TENSORFACES Training Data: In our experiments, we employed gray-level facial training images rendered from 3D scans of 100 subjects. The scans were recorded using a CyberwareTM 3030PS laser scanner and are part of the 3D morphable faces database created at the University of Freiburg [Blanz & Vetter (1999). Each subject was combinatoriall y imaged in Maya from 15 different viewpoints (θ = −60◦ to +60◦ in 10◦ steps on the horizontal plane, φ = 0◦) with 15 different illuminations ( θ = −35◦ to +35◦ in 5◦ increments on a plane inclined at φ = 45◦).
Data Preprocessing: Facial images were warped to an average face template by a piecewise affine transformation given a set of facial landmarks obtained by employing Dlib software [King (2009); Kazemi & Sullivan (2014); Si et al. (2013); Macedo et al. (2006); Hatamizadeh et al. (2019). Illumination was normalized with an adaptive contrast histogram equalization algorithm, but rather than performing contrast correction on the entire image, subtiles of the image were contrast normalized, and tiling artifacts were eliminated through interpolation. Histogram clipping was employed to avoid over-saturated regions.
Experiments: Each image, d ∈ RI0×1, was convolved with five filters banks {Hs‖s = 1...S}. The filtered images, d×0 Hs, resulted in five facial part hierarchies composed of (i) independent pixel parts (ii) parts segmented from different layers of a Gaussian pyramid that were equally or (iii) unequally weighed, (iv) parts were segmented from a Laplacian pyramid that were equally or (v) unequally weighed. We ran five experiments with five facial part hierarchies from which a person representation was computed, Fig. 8. The composite person signature was computed for every test image by employing the multilinear projection algorithm [Vasilescu (2011); Vasilescu & Terzopoulos (2007), and signatures were compared with a nearest neighbor classifier.
To validate the effectiveness of our system on real-world images, we report results on “LFW” dataset (LFW) [Huang et al. (2007). This dataset contains 13,233 facial images of 5,749 people. The photos are unconstrained (i.e., “in the wild”), and include variation due to pose, illumination, expression, and occlusion. The dataset consists of 10 train/test splits of the data. We report the mean accuracy and standard deviation across all splits in Table 2. Fig. 8(b-c) depicts the experimental ROC curves. We follow the supervised “Unrestricted, labeled outside data” paradigm.
Results: While we cannot celebrate closing the gap on human performance, our results are promising. DeepFace, a CNN model, improved the prior art verification rates on LFW from 70% to 97.35%, by training on 4.4M images of 200× 200 pixels from 4, 030 people, the same order of magnitude as the number of people in the LFW database.
We trained on less than one percent (1%) of the 4.4M total images used to train DeepFace. Images were rendered from 3D scans of 100 subjects with an the intraocular distance of approximately 20 pixels and with a facial region captured by 10, 414 pixels (image size ≈ 100× 100 pixels). We have currently achieved verification rates just shy of 80% on LFW. When data is limited, CNN models do not convergence or generalize.
Train on Freiburg: 6 views (±60◦,±30◦,±5◦); 6 illuminations (±60◦,±30◦,±5◦), 45 people Test on Freiburg: 9 views (±50◦, ±40◦, ±20◦, ±10◦, 0◦), 9 illums (±50◦, ±40◦, ±20◦, ±10◦, 0◦), 45 different people LFW Experiment: Models were trained on approximately half of one percent (0.5% < 1%) of the 4.4M images used to train DeepFace. Train on Freiburg: 15 views (±60◦,±50◦, ±40◦,±30◦, ±20◦, ±10◦,±5◦, 0◦), 15 illuminations (±60◦,±50◦, ±40◦,±30◦, ±20◦, ±10◦,±5◦, 0◦), 100 people Test on LFW: We report the mean accuracy and standard deviation across standard literature partitions [Huang et al. (2007), following the Unrestricted, labeled outside data supervised protocol.
Summary: This paper contributes to the tensor algebraic paradigm and models cause-and-effect as a hierarchical block tensor interaction between intrinsic and extrinsic hierarchical causal factors of data formation.
A data tensor expressed as a function of a hierarchical data tensor is a unified tensor model of wholes and parts from which a new compositional hierarchical block tensor factorization was derived. The resulting causal factor representations are interpretable, hierarchical, and statistically invariant to all other causal factors. Our approach was demonstrated in the context of facial images by training on a very small set of synthetic images. While we have not closed the gap on human performance, we report encouraging face verification results on two test data sets–the Freiburg, and the Labeled Faces in the Wild datasets. CNN verification rates improved the 70% prior art to 97.35% when they employed 4.4M images from 4, 030 people, the same order of magnitude as the number of people in the LFW database. We have currently achieved verification rates just shy of 80% on LFW by employing synthetic images from 100 people for a total of less than one percent (1%) of the total images employed by DeepFace. By comparison, when data is limited, CNN models do not convergence, or generalize.

Review summaries:
This paper is an attempt to connect tensor factor analysis with DNN learning. The paper explores the forward causal inference and inverse causal inference problems. The forward causal inference is performed using causal capsules architecture that uses M-mode SVD for subspace learning (autoencoder). The authors show that the M-mode SVD rotates the basis vectors of various clusters such that the causal representation is the same regardless of the data cluster (Fig.1). They introduce causal capsule hierarchy architecture that uses tensor decompositions for forward and inverse causal inference.