Year: 2019

Abstract:
Models that output a vector of responses given some inputs, in the form of a conditional mean vector, are at the core of machine learning. This includes neural networks such as the multilayer perceptron (MLP). However, models that output a symmetric positive definite (SPD) matrix of responses given inputs, in the form of a conditional covariance function, are far less studied, especially within the context of neural networks. Here, we introduce a new variant of the MLP, referred to as the matrix MLP, that is specialized at learning SPD matrices. Our construction not only respects the SPD constraint, but also makes explicit use of it. This translates into a model which effectively performs the task of SPD matrix learning even in scenarios where data are scarce. We present an application of the model in heteroscedastic multivariate regression, including convincing performance on six real-world datasets.

heading: 1 INTRODUCTION
text: For certain applications, it is desirable to construct a conditional covariance matrix as a function of the input (the explanatory variable). The problem arises, for instance, in spatial (and spatio-temporal) statistics, in relation to the heteroscedastic multivariate regression (e.g., Pourahmadi, 1999; Hoff & Niu, 2012; Fox & Dunson, 2015), where we deal with multivariate response measurements for which the typical assumption of homoscedasticity may not be suitable. In such cases, we require models that estimate the covariance matrix that captures the spatial variations in correlations between the elements of the response vector. The covariance matrix is a symmetric positive definite (SPD) matrix which can be challenging to estimate due to its non-Euclidean geometry (Pennec et al., 2005). The central problem that this work is concerned with is learning SPD matrices using neural networks.
To motivate our discussion, consider solving the problem of SPD matrix learning using a multilayer perceptron (MLP) as an example of a fully connected neural network (e.g., Goodfellow et al., 2016). To meet the SPD constraint, one would need to tailor the output layer of the MLP so that the estimated covariance matrix satisfies the SPD requirement. One possible approach would be to use the Cholesky decomposition. The main concern is that this approach does not take into account the non-Euclidean geometry of the SPD matrices. Using empirical evaluations, we will show that the use of “wrong geometry” results in the poor estimation of the SPD matrices in particular where data are scarce.
The primary objective here is to design a nonlinear architecture using neural networks that can effectively perform the task of SPD matrix learning. More specifically, our main contribution is to show how to alter the architecture of the MLP in such a way that it not only respects the SPD constraint, but also makes an explicit use of it. We will achieve this by: 1) Explicitly taking the non-Euclidean geometry of the underlying SPD manifolds (e.g., Pennec et al., 2005) into account by designing a new loss function, and 2) by deriving a new backpropagation algorithm (Rumelhart et al., 1986) that respects the SPD nature of the matrices. This new model will be referred to as the matrix multilayer perceptron (mMLP)1. The mMLP makes use of positive-definite kernels to satisfy the SPD requirement across all layers. Hence, it provides a natural way of enabling deep SPD matrix learning.
We take a step-by-step approach in the development of the model. We first develop a simplified version of the resulting model that is designed for learning SPD matrices (Section 3). We then extend this model into its most general form that can be used for joint estimation of the conditional mean function and the conditional covariance function in a mean-covariance regression setting (Section 3.2). An application of the model is discussed in the context of the heteroscedastic multivariate regression.
1An implementation of this work will be made available via GitHub.

heading: 2 RELATED WORK
text: SPD manifold metric. Earlier approaches for analyzing SPD matrices relied on the Euclidean space. But over the past decade, several studies suggest that non-Euclidean geometries such as the Riemannian structure may be better suited (e.g., Arsigny et al., 2006; Pennec et al., 2005). In this work, we consider the von Neumann divergence (e.g., Nielsen & Chuang, 2000) as our choice of the SPD manifold metric which is related to the Riemannian geometry. Previously, Tsuda et al. (2005) used this divergence in derivation of the matrix exponentiated gradients. Their work suggests its effectiveness for measuring dissimilarities between positive definite (PD) matrices.
SPD manifold learning. There are multiple approaches towards the SPD matrix learning, via the flattening of SPD manifolds through tangent space approximations (e.g., Oncel Tuzel, 2008; Fathy et al., 2016), mapping them into reproducing kernel Hilbert spaces (Harandi et al., 2012; Minh et al., 2014), or geometry-aware SPD matrix learning (Harandi et al., 2014). While these methods typically employ shallow learning, the more recent line of research aims to design a deep architecture to nonlinearly learn target SPD matrices (Ionescu et al., 2015; Huang & Gool, 2017; Masci et al., 2015; Huang et al., 2018). Our method falls in this category but differs in the problem formulation. While the previous methods address the problem where the input is an SPD matrix and the output is a vector, we consider the reverse problem where the input is a matrix with an arbitrary size and the output is an SPD matrix.
Backpropagation. Our extension of the matrix backpropagation differs from the one introduced by Ionescu et al. (2015). In their work, the necessary partial derivatives are computed using a two-step procedure consisting of first computing the functional that describes the variations of the upper layer variables with respect to the variations of the lower layer variables, and then computing the partial derivatives with respect to the lower layer variables using properties of the matrix inner product. In contrast, we make use of the concept of α-derivatives (Magnus, 2010) and its favorable generalization properties to derive a procedure which closely mimics the standard backpropagation.

heading: 3 MATRIX MULTILAYER PERCEPTRON
text: Preliminaries: Matrix α-derivative. Throughout this work we adopt the narrow definition of the matrix derivatives known as the α-derivative (Magnus, 2010) in favor of the broad definition, the ω-derivative. The reason for this is that the α-derivative has better generalization properties. This choice turned out to be crucial in the derivation of the mMLP’s backpropagation routine which involves derivatives of matrix functions w.r.t. the matrix of variables.
Definition: Let F be an m× n matrix function of an n× q matrix of variables X. The α-derivative of F(X) is defined as (Magnus, 2010, Definition 2)
DXF := ∂ vecF(X)
∂ (vecX)> , (1)
where DXF is an mp × nq matrix which contains all the partial derivatives such that each row contains the partial derivatives of one function with respect to all variables, and each column contains the partial derivatives of all functions with respect to one variable.
For convenience, the α-derivative’ basic properties, including the product rule and the chain rule, are summarized in Appendix B.

heading: 3.1 THE BASIC FORM OF THE MMLP
text: Activation matrix function. Let Z = (z1, . . . , zd) denote a matrix of variables zi ∈ Rd. The activation function K(Z) defines a matrix function in the form of [K(Z)]i,j = κ(zi, zj), ∀i, j ∈ {1, . . . , d}, where κ is some differentiable activation function outputting scalar values. In the following, we restrict ourselves to the kernel functions which form PD activation matrix functions. For numerical stability reasons, irrespective of the functional form of κ, we normalize the resulting matrix. This can be achieved by enforcing the trace-one constraint,H(Z) = K(Z)/tr(K(Z)), where H denotes a differentiable SPD activation matrix function of trace one. Without loss of generality, throughout this work, we use the Mercer sigmoid kernel (Carrington et al., 2014) defined as
κ(zi, zj) = tanh(αzi + β) tanh(αzj + β), (2)
where α and β denote the slope and the intercept, respectively. Furthermore, denotes the dot product. In all experiments, we use default values of α = 1 and β = 0. The The α-derivative of the Mercer sigmoid kernel is computed in Appendix E.
Model construction. Let X ∈ Rp1×p2 indicate the input matrix and Y ∈ Rd0×d0 indicate the corresponding output matrix, an SPD matrix of trace one. The mMLP of L hidden layers is shown as mMLP :X→Ŷ and constructed as Ŷ = H(Z0), Z0 = W0H1W>0 + B0, Hl = H(Zl), Zl = WlHl+1W>l + Bl, ∀ 1 ≤ l ≤ L, HL+1 = H(ZL+1), ZL+1 = WL+1vecX(WL+11p1p2)> + BL+1. (3)
The pair of Wl∈Rdl×dl+1 ,∀0 ≤ l≤ L, and WL+1∈RdL+1×p1p2 are the weight matrices, Bl∈Rdl×dl ,∀0≤ l ≤ L+1, are the bias matrices, Zl ∈ Rdl×dl ,∀0 ≤ l≤ L+1, are the latent input matrices, and Hl ∈ Rdl×dl ,∀1 ≤ l ≤ L+ 1, are latent output SPD matrices of trace one. Design choice. In the construction of (3), we have ensured that Hl are SPD matrices of trace one across all layers as opposed to only at the output layer. The idea behind this design choice is to propagate the nonlinearities introduced via the SPD activation matrix functions through all layers. This design choice turned out to be more effective than the alternative, and arguably simpler, design where the SPD requirement is met only at the output layer. We will discuss this further in Section 5.1.3, where we also present an illustrative numerical example.
Loss function. We consider the normalized von Neumann divergence (e.g., Nielsen & Chuang, 2000), also commonly known as the quantum relative entropy (QRE), as the base for the loss function. For two arbitrary SPD matrices of trace one, Φ and Φ̃, the normalized von Neumann divergence is defined as:
∆QRE(Φ̃||Φ) = tr(Φ̃logΦ̃− Φ̃logΦ), (4) where log denotes the matrix logarithm (taking Φ as an example, it is computed using logΦ = Vdiag(logλ)V>, where V and λ are the matrix of eigenvectors and the vector of eigenvalues from the eigendecomposition of Φ). The von Neumann divergence is asymmetric. However, it can be symmetrized by using the fact that the von Neumann entropy of trace one follows the class of generalized quadratic distances (Nielsen & Nock, 2007). Hence, we define the loss function as2
`QRE(Ŷ,Y) = 1
2 (∆QRE(Ŷ||Y) + ∆QRE(Y||Ŷ)). (5)
Taking α-derivative of `QRE involves taking partial derivatives through the eigendecomposition. In the following, we derive a method for analytically computing the derivative of `QRE.
The α-derivative of the symmetrized von Neumann divergence. For the loss function defined in (5), using the α-derivative’s product rule and chain rule (refer to Appendix B), we obtain
DŶ` = 1
2 DŶtr((Ŷ −Y)logŶ)−
1 2 DŶtr(ŶlogY), (6)
where the above two terms are computed using
DŶ(tr((Ŷ −Y)logŶ)) = (vec(logŶ)>)>︸ ︷︷ ︸ 1×d20 +
 vec(Ŷ> −Y>) vec( ∂ ∂Ŷ11 logŶ) vec(Ŷ> −Y>) vec( ∂ ∂Ŷ21 logŶ)
... vec(Ŷ> −Y>) vec( ∂
∂Ŷd0d0 logŶ)
 >
︸ ︷︷ ︸ 1×d20
, (7)
DŶ(tr(ŶlogY)) = (vec(logY) >)>︸ ︷︷ ︸
1×d20
. (8)
2Note that there are multiple ways of symmetrizing the von Neumann divergence. Our choice in (5) resembles to the α-Jensen–Shannon divergence for α = 1 (Nielsen, 2010).
The remaining part in the computation of (6) is to evaluate ∂ ∂Ŷij logŶ, for all i, j ∈ {1, . . . , d0}, which involves taking derivatives through the eigendecomposition. In the following, we take a similar approach as in Papadopoulo & Lourakis (2000) to compute the necessary partial derivatives.
Let Ŷ = Υdiag(λ1, . . . , λd0)Υ > be the eigendecomposition. We can write
∂
∂Ŷij logŶ =
∂
∂Ŷij ΥΛΥ>, where Λ = diag(logλ1, . . . , logλd0),
= ∂Υ
∂Ŷij ΛΥ> + Υ
∂Λ
∂Ŷij Υ> + ΥΛ
∂Υ> ∂Ŷij .
(9)
By multiplying (9) from left and right by Υ> and Υ respectively, we obtain:
Υ> ∂
∂Ŷij logŶ Υ = Υ>
∂Υ
∂Ŷij Λ +
∂Λ
∂Ŷij + Λ
∂Υ> ∂Ŷij Υ = Ξij(Υ)Λ + ∂Λ ∂Ŷij −ΛΞij(Υ), (10)
where we have defined Ξij(Υ) = Υ> ∂∂Ŷij Υ and used the fact that Ξij(Υ) is an antisymmetric matrix, Ξij(Υ) + Ξ>ij(Υ) = 0, which in turn follows from the fact that Υ is an orthonormal matrix,
Υ>Υ = Id0 ⇒ ∂Υ>
∂Ŷij Υ + Υ>
∂Υ
∂Ŷij = Ξ>ij(Υ) + Ξij(Υ) = 0. (11)
Taking the antisymmetric property of Ξij(Υ) into account in (10), we obtain ∂
∂Ŷij logλk = ΥikΥjk, (12)
Ξij(Υkl) = ΥikΥjl + ΥilΥjk 2(logλl − logλk) , ∀l 6= k. (13)
It is notable that by construction, we do not have repeating eigenvalues, that is λk 6= λl, ∀k 6= l, so there exists a unique solution to (13). Once Ξij(Υ) is computed, it follows that
∂Υ
∂Ŷij = ΥΞij(Υ),
∂Υ> ∂Ŷij = −Ξij(Υ)Υ>. (14)
In summary, the necessary partial derivatives for computing (9) is given by (12) and (14). Once (9) is computed for all i, j, we can evaluate (7) and ultimately evalaute (6).
Optimization. The remaining steps are feed-forward computation, backpropagation, and learning, as in the standard MLP. However, here, the backpropagation requires taking derivatives with respect to the matrix functions. These steps are described in Appendix C.

heading: 3.2 THE GENERAL FORM OF THE MMLP
text: We now discuss a general version of the mMLP which produces both a vector and an SPD matrix as outputs. An important application of this model is in heteroscedastic multivariate regression which we will discuss in Section 5.2.
Model construction. As before, let X ∈ Rp1×p2 denote the input matrix. The corresponding outputs in this case are: an SPD matrix of trace one Y ∈ Rd0×d0 and y ∈ Rr0 . The mMLP of L hidden layers is denoted by mMLP :X→{ŷ, Ŷ} and constructed as: ŷ = h(z0), z0 = C0ŶA0h1 + b0, Ŷ = H(Z0), Z0 = W0H1W>0 + B0, hl = h(zl), zl = ClHlAlhl+1 + bl, ∀ 1 ≤ l ≤ L, Hl = H(Zl), Zl = WlHl+1W>l + Bl, ∀ 1 ≤ l ≤ L, hL+1 = h(zL+1), zL+1 = CL+1HL+1AL+11 + bL+1,
HL+1 = H(ZL+1), ZL+1 = WL+1vecX(WL+11p1p2)> + BL+1,
(15)
where hl ∈ Rrl ,Hl ∈ Rdl×dl ,∀1 ≤ l ≤ L+ 1, zl,bl ∈ Rrl , Zl,Bl ∈ Rdl×dl , Cl∈Rrl×dl ,∀0≤ l≤L+1, Al∈Rdl×rl+1, and Wl∈Rdl×dl+1 , ∀0 ≤ l≤ L. Just as in the standard MLP, h is an activation function of choice, e.g., the hyperbolic tangent function.
Loss function. The loss function here needs to be designed with the specific application in mind. In the case of the heteroscedastic multivariate regression, the loss can be defined based on the log-likelihood. This will be discussed in Section 4.
Optimization. The remaining steps of feed-forward computation, backpropagation, and learning, are all described in Appendix D.

heading: 4 THE MMLP IN HETEROSCEDASTIC MULTIVARIATE REGRESSION
text: In this section, we discuss an application of the mMLP in relation to the heteroscedastic multivariate regression, more specifically, the joint mean-covariance regression task.
Model construction. Let Dtrain = {y(i),x(i)}ni=1 be our training dataset consisting of a set of inputs xi ∈ Rdx and a set of responses yi ∈ Rdy . Consider the following multivariate regression problem: y(i) = f(x(i)) + e(i), where f is a nonlinear function, and e(i) is the additive noise on the i-th response measurement. The goal is estimation of the conditional mean function E[y |x∗] and the conditional covariance function Var[y |x∗] for an unseen input x∗ ∈ Dtest. We consider two likelihood models based on our choices of the noise model, namely, the multivariate Gaussian model and its generalization the multivariate power exponential model.
Multivariate Gaussian model. We first consider the noise model to follow a zero-mean multivariate Gaussian (mG) distribution with a dense covariance matrix, that is ei ∼ N (0,Σi). Let Σi = ηiΩi, where ηi = tr(Σi) and Ωi is a trace-one matrix. The noise model can accordingly be reformulated as ei ∼ Ntr1(0,Ωi, ηi) where Ntr1 is referred to as the trace-one mG distribution3. Although these two formulations of the mG distribution are identical, we find it easier to work with the latter. This is because the output layer of the mMLP model in (15) operates under the trace-one constraint. It is of course possible to drop the trace-one constraint from the output layer, but we would then also be forced to use a different kernel function than the Mercer sigmoid kernel in that layer. Instead, we find it easier to work with this reformulated mG distribution with a trace-one covariance matrix which allows us to use the same choice of kernel function (2) across all layers.
Given ei ∼ N (0,Ωi, ηi), the likelihood is defined as `Ntr1 := log Ntr1(y;µ,Ω, η), mMLPθ : x→ {µ,Ω, η}, (16)
where E[y |x] = µ, and Var[y |x] = ηΩ. The set θ includes all the neural network parameters. Multivariate power exponential model. We next consider the noise model to follow a zero-mean multivariate power exponential (mPE) distribution (e.g., Gómez et al., 1998) which is a generalized variant of the mG distribution, that is e(i) ∼ E(0,Σi, αi, βi) where Σi is a symmetric real dispersion matrix and the pair of αi ∈ R+ and βi ∈ R+ control the tail and the shape of the distribution. As a special case, for α = 1 and β = 1, the mPE includes the mG distribution4.
As in the case of the mG distribution, we find it easier to work with a reformulated variant where the dispersion matrix is of trace-one. Let Σi = ηiΩi, where ηi = tr(Σi). The noise model can accordingly be represented as e(i) ∼ Etr1(0,Ωi, αi, βi, ηi) where Etr1 is referred to as the trace-one mPE distribution5.
Given e(i) ∼ Etr1(0,Ωi, αi, βi, ηi), the likelihood is defined as: `Etr1 := log Etr1(y;µ,Ω, α, β, η), mMLPθ : x→ {µ,Ω, α, β, η}, (17)
where E[y |x] = µ and Var[y |x] = ηαν(β)Ω, where ν(β) = 2 1/βΓ( d+12β ) d Γ( d2β ) . The set θ includes all the neural network parameters.
Optimization. Finally, the optimization involves learning the neural network parameters θ by maximizing the likelihoods given by (16) and (17).
3Refer to Appendix F.1 for the exact functional form of the trace-one mG distribution. 4Figure G.1.1 visualizes the probability density function of the distribution for selected values of α and β. 5Refer to Appendix G for the exact functional form of the trace-one mPE distribution and its basic properties.
A. truth D. MLP, B. mMLP, C. mMLP,
Figure 1: SPD matrix learning (refer to Example 1 in Section 5.1.1). Two instances of target covariance (SPD) matrices (20 × 20). (B, C) Estimated covariance matrices by the mMLP using `QRE, and using `quad. (D) Estimated covariance matrices using the Cholesky-based MLP and `quad.

heading: 5 EXPERIMENTS
text: Experiments are divided into two parts. The first part is on the empirical validation of the mMLP model in a supervised task of learning SPD matrices using synthetic data. The second part discusses our results in heteroscedastic multivariate regression on real-world datasets.

heading: 5.1 SPD MATRIX LEARNING

heading: 5.1.1 THE CHOICE OF LOSS FUNCTION
text: Consider the problem of learning SPD matrices on synthetic data using the mMLP model of (3). The objectives are to validate the model, and to evaluate the effect of the choice of the loss function on the overall performance.
The following loss functions are considered for this analysis: The first candidate is the loss function based on the normalized von Neumann divergence `QRE(Ŷ,Y) given by (5). The `QRE is related to the Riemannian geometry. The second candidate is the quadratic loss, `quad(Ŷ,Y) = tr((Ŷ −Y)(Ŷ −Y)>), which is related to the Euclidean geometry.
Example 1. Consider the set Dtrain={xi,Yi}ntraini=1 of inputs xi ∈ R20 and corresponding SPD matrix outputs Yi ∈ Rd0×d0 which are in this case dense covariance matrices (refer to Appendix H.1 for details on the data generation). The goal is to estimate the covariance matrices Ŷ associated to the input vectors from the unseen test set Dtest = {xi}ntesti=1 . The training size is varied between ntrain ={102, 104} samples. The analysis is carried out for d0 ={10, 20}. Two examples of the test target outputs Yi for d0 =20 and ntrain = 102 are visualized in Figure 1-A.
The mMLP models (3) are trained using 3 layers (20 units per layer, dl = 20) under our two choices of loss functions: `QRE, and `quad. All models share the same initialization and the only difference here is the loss function. Refer to Appendix H.1 for additional details on the mMLP initialization 6. The performance is evaluated on the test set, ntest = 103, in terms of the losses as the error measures, shown as EQRE and Equad. Table 1 summarizes the results of the evaluation (also refer to Figure 1-B,C for the visual illustration). The key observation is that the quality of estimates differs considerably depending on the choice of the loss function. The loss function `QRE that takes into account the geometry of the SPD matrices clearly outperforms the one based on the Euclidean geometry, `qaud. The advantage is more pronounced for small set of training data.

heading: 5.1.2 THE MMLP VS. THE MLP BASED ON CHOLESKY DECOMPOSITION
text: As we discussed in Section 1, one can take a heuristic approach and tailor the output layer of a vanilla MLP so that the resulting output matrix satisfies the SPD requirement.
For the sake of comparison, we solve the same task as in Example 1 using the standard MLP with 3 layers (200 units per layer) and with the quadratic loss. To meet the SPD requirement, we use the
6Here and in general throughout this section, special care has been made to minimize the effect of overfitting through trying out various initializations and using early stopping.
Cholesky decomposition at the output layer using the known result that for a SPD matrix there is a unique lower triangular matrix, with ones as its diagonal entires, and a unique diagonal matrix with positive diagonal entires. Table 1 summarizes the results of the evaluation. Overall the performance is quite poor for small set of training data. As the size of training data grows, the performance improves as expected (refer to Figure 1-D for the visual illustration).

heading: 5.1.3 SHALLOW VS DEEP SPD MATRIX LEARNING
text: The design of the mMLP model in (3) enables a mechanism for deep SPD matrix learning by satisfying the SPD constraint across all input, hidden and output layers. The simpler approach would be to consider the standard MLP architecture across input and hidden layers but make use of the activation matrix functions only at the output layer to meet the SPD requirement: Ŷ = H(Z0), Z0 = W0h1(W01)> + B0, hl = h(zl), zl = Wlhl+1 + bl, 1 ≤ l ≤ L, hL+1 = h(zL+1), zL+1 = WL+1vecX + bL+1. (18)
This amounts to a shallow design in the sense that it does not enable a mechanism for preserving the SPD constraint across all layers during the learning.
The design in (3) allows nonlinearities to pass through layers via activation function matrices imposing the SPD constraint, whereas in the shallow design, nonlinearities are propagated across layers via activation functions without imposing any constraints. Our hypothesis is that the former has advantage over the latter in that it captures complex dependencies which are important for the SPD matrix learning. Below, we present a numerical example which indeed highlights the importance of preserving the SPD constraint across all layers when learning the SPD matrix.
Example 2. Consider a similar experiment as in Example 1 for the case of ntrain = 102 and output dimensions d0 = {10, 20}. We directly compare the performance of (3) against (18) under different number of hidden layers L = {2, 4, 6}. For the mMLP model, the number of hidden units at each layer is set to 20, and for the MLP model it is set to 200 units. The shallow design (18) uses the hyperbolic tangent as the activation function h(·). The same choice of the activation matrix function H(·), given by (2), is used for both models. We use `QRE as the choice of the loss function for both models (refer to Appendix H.2 for additional details on the initialization). The performance is evaluated in terms of EQRE.
Table 2 summarizes the results of the evaluation. Although the shallow design (18) performs relatively well, it underperforms in comparison to (3). Given the limited number of training samples, arbitrarily increasing the number of layers may not be necessarily advantageous, which is the case for both models. However, in this regard, the design in (18) is more sensitive.

heading: 5.2 EXPERIMENTAL RESULTS ON THE MULTI-OUTPUT REGRESSION DATASETS
text: In heteroscedastic multivariate regression, the central task is to estimate the conditional covariance matrix that captures the spatial variations in correlations between the elements of the response vector. The underlying hypothesis is that if a model can effectively capture the heteroscedasticity, then it will provide better uncertainty quantification, and the estimation of the conditional mean response should also improve, in comparison to the model that is build based on the homoscedasticity assumption.
For this purpose, we compare our mMLP-based regression models against another heteroscedastibased mean-covariance regression model by Fox & Dunson (2015) and two homoscedastic-based mean-regression models, namely, the MLP regressor and the Gaussian process (GP). The models used in this experiment are summarized in Table 3. Real-world datasets. We compare the performance of the heteroscedastic and homoscedastic regression models on six real-world multi-output datasets. Key features of the datasets are summarized in Table 5a. The performance is evaluated in terms of the root-mean-square error (RMSE) on test sets, shown in Table 5b. The result suggests that the mMLP homoscedastic-based regression models are capable of capturing dependencies between the output measurements which contributes to the better estimation of the mean predictions.

heading: 6 LIMITATIONS AND FUTURE WORK
text: The main limitation of the mMLP has to do with scalability to higher dimensions. The complexity associated with computing the α-derivative of the von Neumann loss function (5) at the output layer is O(d30). Taking the symmetric nature of the SPD matrices into account, the computational complexity at the hidden layer l reduces to O(d2l ). Our implementation of the matrix backpropagation involves the use of multiple Kronecker products. Although it facilitates the implementation, we would need access to the full Jacobian matrices (d2l × d2l ). However, these matrices are in fact available in the form of sparse block matrices, which means that it is possible to implement a memory-efficient computation of the tensor products without the need to actually have access to the full matrices. Future work is needed in this direction.
We believe that there are many cases in which the mMLP can prove to be useful. An interesting direction for future work is to investigate application of the model in the context of the conditional density estimation within the framework of variational autoencoders.

heading: 7 DISCUSSION
text: We introduced a new method to learn SPD matrices, referred to as the matrix multilayer perceptron (mMLP). The mMLP takes the non-Euclidean geometry of the underlying SPD manifolds into account by making use of the von Neumann divergence as the choice of the SPD manifold metric. One key aspect of the mMLP is that it preserves the SPD constraint across all layers by exploiting PD kernel functions and a backpropagation algorithm that respects the inherent SPD nature of the matrices. We studied application of the model in the context of heteroscedastic multivariate regression. We showed the effectiveness of the proposed model on multiple real-world datasets.

heading: A MATRIX NOTATIONS
text: We use > for the transpose operator, tr(·) for the trace operator, and det(·) for the matrix determinant. The symmetric part of a square matrix B is denoted by sym(B) = (B + B>)/2. The Kronecker product is denoted by ⊗, the Hadamard product by ◦, and the dot product by . We use the vec-operator for column-by-column stacking of a matrix A, shown as vecA ≡ vec(A). Let A be an m× n matrix, the operator P(m,n) will then rearrange vecA to its matrix form: A = P(m,n)(vecA). For the m× n dimensional matrix A, the commutation matrix is shown asK(m,n) which is the mn×mn matrix that transforms vecA into vecA> as: K(m,n)vecA = vecA>. An m×m identity matrix is shown as Im. If T (X) : Rd×d → R is a real-valued function on matrices, then ∇XT (X) denotes the gradient with respect to the matrix X,∇XT (X) = [ ∂T ∂xij ] i,j=1:d
. The matrix logarithm and the matrix exponential are written as logA and expA, respectively. The matrix exponential in the case of symmetric matrices can be expressed using the eigenvalue decomposition as expA = V(expΛ)V>, where V is an orthonormal matrix of eigenvectors and Λ is a diagonal matrix with the eigenvalues on the diagonal. The matrix logarithm is the inverse of the matrix exponential if it exists. If A is symmetric and strictly positive definite (PD) it is computed using logA = V(logΛ)V>, where (logΛ)i,i = logΛi,i.
B THE α-DERIVATIVE: DEFINITION AND PROPERTIES
Definition Let F be an m× n matrix function of an n× q matrix of variables X. The α-derivative of F(X) is defined as (Magnus, 2010, Definition 2)
DXF := ∂ vecF(X)
∂ (vecX)> , (19)
where DXF is an mp × nq matrix which contains all the partial derivatives such that each row contains the partial derivatives of one function with respect to all variables, and each column contains the partial derivatives of all functions with respect to one variable.
Product rule Let F (m× p) and G (p× r) be functions of X (n× q). Then the product rule for the α-derivative is given by (Magnus, 2010)
DX(FG) = (G > ⊗ Im)DXF + (Ir ⊗ F)DXG. (20)
Chain rule Let F (m× p) be differentiable at X (n× q), and G (l × r) be differentiable at Y = F(X), then the composite function H(X) = G(F(X)) is differentiable at X, and
DXH = DYGDXF, (21)
which expresses the chain rule for the α-derivative (Magnus, 2010).

heading: C THE BASIC CASE OF THE MMLP

heading: C.1 FEEDFORWARD STEP
text: At the feed-forward computation, we compute and store the latent outputs Ŷ, Hl for all l ∈ {L+ 1, . . . , 1} using the current setting of the parameters, which are Wl, and Bl computed from the learning step, Appendix C.3.

heading: C.2 BACKPROPAGATION STEP
text: We first summarize the necessary α-derivatives for the backpropagation, and then write down the backpropagation procedure accordingly.
Derivatives required for backpropagation The derivative of the activation matrix function depends on the specific choice of kernel function, and in general it is computed readily from the
definition of α-derivative,
DZlH(Zl) := ∂ vecH(Zl) ∂ (vecZl)> , l ∈ {0, . . . , L+ 1}. (22)
For our specific choice of activation function, the Mercer Sigmoid kernel (2), it is computed in Appendix E.
Via repeated use of the product rule of α-derivatives (20), we obtain
DWlZl = (Wl ⊗ Idl)(H>l+1 ⊗ Idl) + (Idl ⊗ (WlHl+1))K(dl,dl+1), l ∈ {0, . . . , L}, (23) DWL+1Zl = (WL+11p1p2 ⊗ IdL+1)((vecX)> ⊗ IdL+1) + (IdL+1 ⊗ (WL+1vecX))(IdL+1 ⊗ 1>p1p2)K(dL+1,p1p2), (24) DHl+1Zl = (Wl ⊗ Idl)(Idl+1 ⊗Wl), l ∈ {0, . . . , L}, (25)
whereK is the commutation matrix (refer to Appendix A for a summary of the matrix notation).
Backpropagation In the interest of simple expressions, let H0 ≡ Ŷ. Backpropagation to the hidden layer l is computed recursively using the derivatives computed at the previous layer according to
DZl` = DHl`DZlHl, ∀l ∈ {0, . . . , L+ 1}, (26) DWl` = DZl`DWlZl, ∀l ∈ {0, . . . , L+ 1}, (27) DHl+1` = DZl`DHl+1Zl, ∀l ∈ {0, . . . , L}, (28) DBl` = DZl`, ∀l ∈ {0, . . . , L+ 1}. (29)

heading: C.3 LEARNING STEP
text: Learning involves updating the weights Wl and the biases Bl using derivatives computed during the backpropagation step. These are updated using derivatives DWl` and DBl` for a given learning rate η as
Wl ←Wl − ηP(dl,dl+1)(DWl`), ∀l ∈ {0, . . . , L+ 1}, (30) Bl ← Bl − ηP(dl,dl)(DBl`), ∀l ∈ {0, . . . , L+ 1}, (31)
where P is the rearrangement operator introduced in Appendix A.

heading: D THE GENERAL FORM OF THE MMLP

heading: D.1 FEEDFORWARD STEP
text: The forward path involves computing and storing both hl, ŷ and Hl, Ŷ using the current settings of the parameters, for all l ∈ {0, . . . , L+ 1}.

heading: D.2 BACKPROPAGATION STEP
text: Most of the necessary derivatives are identical to the ones in Appendix C.2. However, there are some additional derivatives needed which we will discuss in the following. We then write down the backpropagation formula.

heading: D.2.1 REQUIRED DERIVATIVES FOR BACKPROPAGATION
text: The derivative of the activation function depends on the choice of the function, and it is computed using the definition of the α-derivative,
Dzlhl(zl) = ∂ h(zl)
∂ (zl)> , l ∈ {0, . . . , L+ 1}. (32)
The other required derivatives are computed as
DAlzl = ClHl, l ∈ {0, . . . , L}, (33) DAL+1zL+1 = (CL+1HL+1)(1 > rL+1 ⊗ IdL+1), (34) DClzl = (HlAlhl+1) > ⊗ Irl , l ∈ {0, . . . , L}, (35) DCL+1zL+1 = (HL+1AL+11rL+1) > ⊗ IrL+1 , (36) Dhl+1zl = ClHlAl, l ∈ {0, . . . , L}, (37) DHlzl = ((Alhl+1)
> ⊗ Irl)(Idl ⊗Cl), l ∈ {0, . . . , L}, (38) DHL+1zL+1 = ((AL+11L+1) > ⊗ IrL+1)(IdL+1 ⊗CL+1). (39)
Backpropagation For simplicity of expressions, let h0 ≡ ŷ and H0 ≡ Ŷ. The derivatives are recursively computed as
Dh0` ≡ Dŷ` (40) DH0` ≡ DŶ` = Dz0`DŶz0 + DŶ` (41) Dzl` = Dhl`Dzlhl, ∀l ∈ {0, . . . , L+ 1}, (42) Dhl+1` = Dzl`Dhl+1zl, ∀l ∈ {0, . . . , L}, (43) DZl` = DHl`DZlHl, ∀l ∈ {0, . . . , L+ 1}, (44) DHl+1` = Dzl+1`DHl+1zl+1 + DZl`DHl+1Zl, ∀l ∈ {0, . . . , L}, (45) DAl` = Dzl`DAlzl, ∀l ∈ {0, . . . , L+ 1}, (46) DCl` = Dzl`DClzl, ∀l ∈ {0, . . . , L+ 1}, (47) DWl` = DZl`DWlZl, ∀l ∈ {0, . . . , L+ 1}, (48) Dbl` = Dzl`, ∀l ∈ {0, . . . , L+ 1}, (49) DBl` = DZl`, ∀l ∈ {0, . . . , L+ 1}. (50)

heading: D.3 LEARNING STEP
text: The learning step involves updating the weights and the biases which are computed using derivatives computed from the backpropagation step. Update rules for Wl, and Bl are the same as the ones given in Appendix C.3. The remaining parameters are learned in a similar fashion,
Al ← Al − ηP(dl,rl+1)(DAl`), ∀l ∈ {0, . . . , L+ 1}, (51) Cl ← Cl − ηP(rl,dl)(DCl`), ∀l ∈ {0, . . . , L+ 1}, (52) bl ← bl − ηP(rl,1)(Dbl`), ∀l ∈ {0, . . . , L+ 1}. (53)
E THE α-DERIVATIVE OF THE MERCER SIGMOID KERNEL
The α-derivative of the Mercer sigmoid kernel can be computed as
DZH =  · · · · · · · · · · · · ∂ ∂zi κmn trK︸ ︷︷ ︸
1×dl
· · ·
· · · · · · · · ·  , ∀ i,m, n ∈ {1, . . . , dl}, (54) where zi indicates the ith column of Z, trK ≡ trK(Z), κmn ≡ κ(zm, zn) as defined in (2), and
∂
∂zi
κmn trK =
( α(1> − f(zi) ◦ f(zi)) (trK)2 ) ◦ ( trK f ( ∂ ∂zi (zm ◦ zn) ) − 2κmnf(zi) ) , (55)
where f(zi) := tanh(αzi − β).
Under review as a conference paper at ICLR 2020
13 Supplementary figures195
Figure 1: The probability density function of a trace-one mPE distribution, Etr1(µ, M/⌘, ⌘,↵, ) for fixed values of µ, M , ⌘ = tr(M) and varying values of scale and shape parameters (↵, ). When ↵ = 1 and = 1, the density corresponds to the multivariate Gaussian distribution N (µ, M).
10
13 Supplementary figures195
Figure 1: The probability density function of a trace-one mPE distribution, Etr1(µ, M/⌘, ⌘,↵, ) for fixed values of µ, M , ⌘ = tr(M) and varying values of scale and shape parameters (↵, ). When ↵ = 1 and = 1, the density corresponds to the multivariate Gaussian distribution N (µ, M).
10
13 Supplementary figures195
Figure 1: The probability density function of a trace-one mPE distribution, Etr1(µ, M/⌘, ⌘,↵, ) for fixed values of µ, M , ⌘ = tr(M) and varying values of scale and shape parameters (↵, ). When ↵ = 1 and = 1, the density corresponds to the multivariate Gaussian distribution N (µ, M).
10
α
β
Figure G.1: The probability density function of a trace-one mPE distribution, Etr1(µ,M/η, η, α, β) for fixed values of µ,M , η = tr(M) and varying values of scale and shape parameters (α, β). When α = 1 and β = 1, the density corresponds to the multivariate Gaussian distribution N (µ,M).

heading: F TRACE-ONE MULTIVARIATE GAUSSIAN DISTRIBUTION

heading: F.1 PROBABILITY DENSITY FUNCTION
text: For a d-dimensional random variable ϑ ∈ Rd, we define the trace-one Gaussian distribution according to
Ntr1(ϑ;µ,Ω,η)= 1
det(2πηΩ) 1 2
e− 1 2 (ϑ−µ) >(ηΩ)−1(ϑ−µ), (56)
where µ ∈ Rd is the mean, η ∈ R+ is the scale parameter, and Ω ∈ Rd×d is the trace-one covariance matrix, tr(Ω) = 1.
F.2 THE α-DERIVATIVES OF THE LOG-PDF
The α-derivatives of the trace-one Gaussian distribution’s log-pdf with respect to its parameters are summarized as
DΩlogNtr1(ϑ;µ,Ω, η) = − 1 2 (vec(Ω−1))> − 1 2 (DΩt) >, (57) DΩt = −vec((ηΩ)−1(ϑ− µ)(ϑ− µ)>Ω−1), (58)
DµlogNtr1(ϑ;µ,Ω, η) = (ϑ− µ)(ηΩ)−1, (59)
DlogηlogNtr1(ϑ;µ,Ω, η) = − d
2 +
1 2 (ϑ− µ)>(ηΩ)−1(ϑ− µ). (60)

heading: G TRACE-ONE MULTIVARIATE POWER EXPONENTIAL (MPE) DISTRIBUTION

heading: G.1 PROBABILITY DENSITY FUNCTION
text: For a random variable ϑ ∈ Rd, the trace-one mPE distribution can be expressed as
Etr1(ϑ;µ,Ω, α, β, η)= c(α, β)
(det(ηΩ)) 1 2
exp { −1
2
( t(ϑ;µ,Ω)
αη
)β} ,
c(α, β) = βΓ(d2 )
π d 2 Γ( d2β )2 d 2β α d 2
, t(ϑ;µ,Ω) := (ϑ− µ)>Ω−1(ϑ− µ), tr(Ω) = 1. (61)
Here, µ is the mean vector, and Ω is a d× d symmetric real dispersion matrix where tr(Ω) = 1. The parameter η has the same role as in the trace-one Gaussian distribution. The pair of α ∈ R+
and β ∈ R+ control the tail and the shape of the distribution. As a special case, the mPE includes the Gaussian distribution: For α = 1 and β = 1, the trace-one mPE distribution corresponds to the trace-one multivariate Gaussian distribution. Figure G.1.1 visualizes the probability density function of the distribution for selected values of α and β, for the case of d = 2.
Remark: Very large and very small values of α and β might be undesirable as, for one, they pose numerical challenges. In practice, these parameters can be restricted within a range by choosing suitable output activation functions. In all experiments in this paper, we choose to bound them conservatively as: 0.5 ≤ α, β ≤ 1.5, using the sigmoid function.

heading: G.2 MOMENTS
text: Let ϑ ∈ Rd and ϑ ∼ Etr1(µ,Ω, α, β, η). The mPE’s mean vector and covariance matrix are computed from:
E[ϑ] = µ, V[ϑ] = αην(β)Ω, ν(β) := 21/βΓ(d+22β )
dΓ( d2β ) , (62)
where Γ(·) denotes the gamma function.
G.3 THE α-DERIVATIVES OF THE LOG-PDF
It is straightforward to take derivatives of the mPE’s log-pdf using the favorable generalization properties of the α-derivative’s chain and product rules. These are summarized as:
DΩlogEtr1(ϑ;µ,Ω, α, β, η) = − 1 2 (vec(Ω−>))> − β 2αη (t/αη)β−1DΩt, (63)
DΩt = −(vec(Ω−1(ϑ− µ)(ϑ− µ)>Ω−1)>)>, (64)
DµlogEtr1(ϑ;µ,Ω, α, β, η) = − β
2αη (t/αη)β−1Dµt, (65)
Dµt = −2(vec((ϑ− µ)>Ω−1)>)>, (66)
DlogηlogEtr1(ϑ;µ,Ω, α, β, η) = − d
2 +
βt
2αη (t/αη)β−1, (67)
DαlogEtr1(ϑ;µ,Ω, α, β, η) = − d
2α +
βt
2ηα2 (t/αη)β−1, (68)
Dβ logEtr1(ϑ;µ,Ω, α, β, η) = Dβ logc(α, β)− 1
2 (t/αη)β log(t/αη), (69)
Dβ logc(α, β) = 1
β +
d
2β2 (ψ(d/2β) + log2). (70)

heading: H ADDITIONAL DETAILS ON THE EXPERIMENTS

heading: H.1 EXAMPLE 1
text: Data generation Let A ∈ Rd0×d0 be a matrix where each of its elements is generated from a standard normal distribution. The matrixA is kept fixed. The ith class covariance Yi is computed according to the following procedure:
1. Draw 104 samples from a known Gaussian distribution N (µi,Σi) with a unique mean µi ∈ Rd0 and a unique dense covariance matrix Σi ∈ Rd0×d0 .
2. Let tj be a random sample from this Gaussian. For this sample, compute yj = Atj . For all 104 samples, collect the results into y = {yj}10 4
j=1. 3. Compute the sample covariance of y and normalize the resulting covariance matrix to trace
one, that is Yi ← cov(y)/tr(cov(y)).
Initialization All models use the same batch size (equal to 5), the same choice of activation matrix function, which is given by the Mercer sigmoid kernel (2), and the same optimizer (the Adam optimizer (Kingma & Welling, 2014) with default settings).

heading: H.2 EXAMPLE 2
text: Data generation See the data generation procedure in Example 1.
Initialization Both models (18) and (3) use the same batch size (equal to 5), the same choice of loss function (5), and the same optimizer (the Adam optimizer (Kingma & Welling, 2014) with default settings). Both models use the same choice of the output activation matrix function, given by the Mercer sigmoid kernel (2). The model in (18) uses the hyperbolic tangent as the activation function across the hidden layers, while (3) makes use of the same choice of the activation matrix function as in its output layer.

heading: H.3 MULTI-OUTPUT REGRESSION DATASETS
text: H.3.1 oes10
The dataset oes10 was obtained from Spyromitros-Xioufis et al. (2016). The Occupational Employment Survey (OES) datasets contain records from the years of 2010 (OES10) of the annual Occupational Employment Survey compiled by the US Bureau of Labor Statistics. As described in (Spyromitros-Xioufis et al., 2016), "each row provides the estimated number of full-time equivalent employees across many employment types for a specific metropolitan area". We selected the same 16 target variables as listed in (Spyromitros-Xioufis et al., 2016, Table 5). The remaining 298 variables serve as the inputs. Data samples were randomly divided into training and test sets (refer to Table 5a).
H.3.2 edm
The dataset edm was obtained from Karalic & Bratko (1997). The electrical discharge machining (EDM) dataset contain domain measurements in which the workpiece surface is machined by electrical discharges occurring in the gap between two electrodes: the tool and the workpiece. Given the two input variables, gap control and flow control, the aim here is to predict the other 16 target variables representing mean values and deviations of the observed quantities of the considered machining parameters. Data samples were randomly divided into training and test sets (refer to Table 5a).
H.3.3 atp1d AND atp7d
The datasets atp1d and atp7d were obtained from (Spyromitros-Xioufis et al., 2016). The Airline Ticket Price (ATP) dataset includes the prediction of airline ticket prices. As described in (Spyromitros-Xioufis et al., 2016), the target variables are either the next day price, atp1d, or minimum price observed over the next seven days atp7d for 6 target flight preferences listed in (Spyromitros-Xioufis et al., 2016, Table 5). There are 411 input variables in each case. The inputs for each sample are values considered to be useful for prediction of the airline ticket prices for a specific departure date, for example, the number of days between the observation date and the departure date, or the boolean variables for day-of-the-week of the observation date. Data samples were randomly divided into training and test sets (refer to Table 5a).
H.3.4 scm1d AND scm20d
The datasets scm1d and scm20d were obtained from (Spyromitros-Xioufis et al., 2016). The Supply Chain Management (SCM) datasets are derived from the Trading Agent Competition in Supply Chain Management (TAC SCM) tournament from 2010. As described in (Spyromitros-Xioufis et al., 2016), each row corresponds to an observation day in the tournament. There are 280 input variables in these datasets which are observed prices for a specific tournament day. The datasets contain 16 regression targets, where each target corresponds to the next day mean price scm1d or mean price for 20 days in the future scm20d for each product (Spyromitros-Xioufis et al., 2016, Table 5). Data samples were randomly divided into training and test sets (refer to Table 5a).

Review summaries:
This manuscript proposes a novel formulation of the MLP to address predicting a symmetric positive definite (SPD) matrix from an input vector or matrix.  While the field has had methods for years to estimate SPD matrices (such as the covariance matrix estimate in the reparameterization trick), this manuscript proposes a markedly different approach based on a different layer structure and repeated normalization steps based on Mercer Kernels.  Additionally, the loss used can be modified to use more PSD-specific losses, such as the symmetrized von Neumann divergence, rather than the traditional quadratic loss.  This loss appears to give significantly better solutions on synthetic data.

This paper explores the problem of deep heteroskedastic multivariate regression where the goal is to regress over symmetric positive definite matrices; that is, the deep learning model should take as input data points, and produce a conditional covariance matrix as the output. The key challenge in this setting is how to ensure the predicted matrix is positive definite (and thus follows the non-linear geometry of these matrices), how the neural network can be trained for this task, and what loss function can be used for effective training. The paper proposes a neural network with bilinear layers in this regard, and uses the von Neumann divergence as the loss function to regress the predicted covariance against a ground truth SPD matrix. The gradients of the von Neumann divergence are provided for learning via backpropagation. Experiments on several synthetic datasets and small scale datasets are provided, showcasing some benefits.