Title: MATRIX MULTILAYER PERCEPTRON

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.

1. INTRODUCTION:

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.

2. RELATED WORK:

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.

3. MATRIX MULTILAYER PERCEPTRON:

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)

EQUATION

where D X F 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.

3.1. THE BASIC FORM OF THE MMLP:

Activation matrix function.

Let Z = (z 1 , . . . , z d ) denote a matrix of variables z i ∈ R d .

The activation function K(Z) defines a matrix function in the form of [K(Z)] i,j = κ(z i , z j ), ∀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

EQUATION

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 ∈ R p1×p2 indicate the input matrix and Y ∈ R d0×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(Z 0 ), Z 0 = W 0 H 1 W 0 + B 0 , H l = H(Z l ), Z l = W l H l+1 W l + B l , ∀ 1 ≤ l ≤ L, H L+1 = H(Z L+1 ), Z L+1 = W L+1 vecX(W L+1 1 p1p2 ) + B L+1 .

(3)

The pair of W l ∈ R d l ×d l+1 , ∀0 ≤ l ≤ L, and W L+1 ∈ R d L+1 ×p1p2 are the weight matrices, B l ∈ R d l ×d l , ∀0 ≤ l ≤ L+1, are the bias matrices, Z l ∈ R d l ×d l , ∀0 ≤ l ≤ L+1, are the latent input matrices, and H l ∈ R d l ×d l , ∀1 ≤ l ≤ L + 1, are latent output SPD matrices of trace one.

Design choice. In the construction of (3), we have ensured that H l 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:

EQUATION

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 

EQUATION

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

EQUATION

where the above two terms are computed using

EQUATION

)

The remaining part in the computation of ( 6) is to evaluate ∂ ∂ Yij log Y, for all i, j ∈ {1, . . . , d 0 }, 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

EQUATION

By multiplying ( 9) from left and right by Υ and Υ respectively, we obtain:

EQUATION

where we have defined

Ξ ij (Υ) = Υ ∂ ∂ Yij Υ 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,

EQUATION

)

Taking the antisymmetric property of Ξ ij (Υ) into account in (10), we obtain

EQUATION

EQUATION

It is notable that by construction, we do not have repeating eigenvalues, that is λ k = λ l , ∀k = l, so there exists a unique solution to (13

EQUATION

)

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.

3.2. THE GENERAL FORM OF THE MMLP:

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 ∈ R p1×p2 denote the input matrix. The corresponding outputs in this case are: an SPD matrix of trace one Y ∈ R d0×d0 and y ∈ R r0 . The mMLP of L hidden layers is denoted by mMLP : X → { y, Y} and constructed as:

EQUATION

where

h l ∈ R r l , H l ∈ R d l ×d l , ∀1 ≤ l ≤ L + 1, z l , b l ∈ R r l , Z l , B l ∈ R d l ×d l , C l ∈ R r l ×d l ,∀0 ≤ l ≤ L+1, A l ∈ R d l ×r l+1 , and W l ∈ R d l ×d l+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.

4. THE MMLP IN HETEROSCEDASTIC MULTIVARIATE REGRESSION:

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 D train = {y (i) , x (i) } n i=1 be our training dataset consisting of a set of inputs x i ∈ R dx and a set of responses y i ∈ R dy . 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 * ∈ D test .

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 e i ∼ 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 e i ∼ N tr1 (0, Ω i , η i ) where N tr1 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 e i ∼ N (0, Ω i , η i ), the likelihood is defined as

EQUATION

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) ∼ E tr1 (0, Ω i , α i , β i , η i ) where E tr1 is referred to as the trace-one mPE distribution5 . Given e (i) ∼ E tr1 (0, Ω i , α i , β i , η i ), the likelihood is defined as:

EQUATION

where

E[y | x] = µ and Var[y | x] = ηαν(β)Ω, where ν(β) = 2 1/β Γ( d+1 2β ) d Γ( d 2β )

. 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). The mMLP models (3) are trained using 3 layers (20 units per layer, d l = 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 initialization6 . The performance is evaluated on the test set, n test = 10 3 , in terms of the losses as the error measures, shown as E QRE and E quad . 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.

5.1.2. THE MMLP VS. THE MLP BASED ON CHOLESKY DECOMPOSITION:

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 

d 0 = 20, n train = 10 4 E quad E QRE 3.3×10 -8 7.7×10 -5 2.1 × 10 -3 3.9 × 10 -2 1.3 × 10 -2 8.1 × 10 -1 d 0 = 10, n train = 10 2 E quad E QRE

6.3×10 -6 1.1×10 -4 3.7 × 10 -2 6.6 × 10 -1 6.5 × 10 -1 4.5

d 0 = 20, n train = 10 2 E quad E QRE 4.3×10 -5 2.2×10 -3 6.0 × 10 -2 1.2 7.3 × 10 -1 6.4

Table 2 : Shallow vs deep architecture, Example 2. Error is measured in terms of E QRE .

d 0 = 10 d 0 = 20 Model L = 2 L = 4 L = 6 L = 2 L = 4 L = 6 Deep architecture (3) 2×10 -4 1×10 -5 4×10 -5 5×10 -3 3×10 -4 6 × 10 -4 Shallow architecture (18) 4×10 -3 5×10 -3 4×10 -2 8×10 -2 6×10 -2 7×10 -1

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).

5.1.3. SHALLOW VS DEEP SPD MATRIX LEARNING:

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:

EQUATION

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 n train = 10 2 and output dimensions d 0 = {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 E QRE .

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.

5.2. EXPERIMENTAL RESULTS ON THE MULTI-OUTPUT REGRESSION DATASETS:

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. 

H(Z) = K(Z) tr(K(Z))

, where [K(Z)] i,j = κ(z i , z j ), ∀i, j. The kernel function κ(•, •) is given by (2). 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.

6. LIMITATIONS AND FUTURE WORK:

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(d 3 0 ). Taking the symmetric nature of the SPD matrices into account, the computational complexity at the hidden layer l reduces to O(d 2 l ). 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 (d 2 l × d 2 l ). 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.

7. DISCUSSION:

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.

A MATRIX NOTATIONS:

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 as K (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 I m . If T (X) : R d×d → R is a real-valued function on matrices, then ∇ X T (X) denotes the gradient with respect to the matrix X, ∇ X T (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)

EQUATION

where D X F 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 )

EQUATION

)

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

EQUATION

which expresses the chain rule for the α-derivative (Magnus, 2010) .

C THE BASIC CASE OF THE MMLP:

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

C.2 BACKPROPAGATION STEP:

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,

EQUATION

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

EQUATION

EQUATION

EQUATION

where K is the commutation matrix (refer to Appendix A for a summary of the matrix notation).

Backpropagation In the interest of simple expressions, let H 0 ≡ Y. Backpropagation to the hidden layer l is computed recursively using the derivatives computed at the previous layer according to

EQUATION

EQUATION

EQUATION

EQUATION

C.3 LEARNING STEP Learning involves updating the weights W l and the biases B l using derivatives computed during the backpropagation step. These are updated using derivatives D W l and D B l for a given learning rate η as

EQUATION

EQUATION

)

where P is the rearrangement operator introduced in Appendix A.

D THE GENERAL FORM OF THE MMLP D.1 FEEDFORWARD STEP:

The forward path involves computing and storing both h l , y and H l , Y using the current settings of the parameters, for all l ∈ {0, . . . , L + 1}.

D.2 BACKPROPAGATION STEP:

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.

D.2.1 REQUIRED DERIVATIVES FOR BACKPROPAGATION:

The derivative of the activation function depends on the choice of the function, and it is computed using the definition of the α-derivative,

EQUATION

The other required derivatives are computed as

EQUATION

EQUATION

EQUATION

EQUATION

EQUATION

EQUATION

EQUATION

Backpropagation For simplicity of expressions, let h 0 ≡ y and H 0 ≡ Y. The derivatives are recursively computed as

EQUATION

EQUATION

EQUATION

EQUATION

EQUATION

EQUATION

EQUATION

EQUATION

EQUATION

D.3 LEARNING STEP:

The learning step involves updating the weights and the biases which are computed using derivatives computed from the backpropagation step. Update rules for W l , and B l are the same as the ones given in Appendix C.3. The remaining parameters are learned in a similar fashion,

EQUATION

EQUATION

EQUATION

E THE α-DERIVATIVE OF THE MERCER SIGMOID KERNEL:

The α-derivative of the Mercer sigmoid kernel can be computed as

EQUATION

where z i indicates the i th column of Z, trK ≡ trK(Z), κ mn ≡ κ(z m , z n ) as defined in (2), and

EQUATION

where f(z i ) := tanh(αz iβ). 

EQUATION

where µ ∈ R d is the mean, η ∈ R + is the scale parameter, and Ω ∈ R d×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

EQUATION

D Ω t = -vec((ηΩ) -1 (ϑµ)(ϑµ) Ω -1 ), 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.

G.2 MOMENTS:

Let ϑ ∈ R d and ϑ ∼ E tr1 (µ, Ω, α, β, η). The mPE's mean vector and covariance matrix are computed from:

EQUATION

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: Data generation Let A ∈ R d0×d0 be a matrix where each of its elements is generated from a standard normal distribution. The matrix A is kept fixed. The i th class covariance Y i is computed according to the following procedure:

EQUATION

EQUATION

1. Draw 10 4 samples from a known Gaussian distribution N (µ i , Σ i ) with a unique mean µ i ∈ R d0 and a unique dense covariance matrix Σ i ∈ R d0×d0 . 2. Let t j be a random sample from this Gaussian. For this sample, compute y j = At j . For all 10 4 samples, collect the results into y = {y j } 10 4 j=1 . 3. Compute the sample covariance of y and normalize the resulting covariance matrix to trace one, that is Y i ← cov(y)/tr(cov(y)).

:

An implementation of this work will be made available via GitHub.

Note 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).

Refer to Appendix F.1 for the exact functional form of the trace-one mG distribution.

Figure G.1.1 visualizes the probability density function of the distribution for selected values of α and β.

Refer to Appendix G for the exact functional form of the trace-one mPE distribution and its basic properties.

Here 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.

Figures and Tables:

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 .

Figure 1: The probability density function of a trace-one mPE distribution, E tr1 (µ, 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 ).

Figure 1: The probability density function of a trace-one mPE distribution, E tr1 (µ, 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 ).

µ logN tr1 (ϑ; µ, Ω, η) = (ϑµ)(ηΩ) -1 , (59)D logη logN tr1 (ϑ; µ, Ω, η) =variable ϑ ∈ R d, the trace-one mPE distribution can be expressed asE tr1 (ϑ; µ, Ω, α, β, ηϑ; µ, Ω) := (ϑµ) Ω -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.

µ logE tr1 (ϑ; µ, Ω, α, β, η) = -β 2αη (t/αη) β-1 D µ t,(65)D µ t = -2(vec((ϑµ) Ω -1 ) ) ,(66)D logη logE tr1 (ϑ; µ, Ω, α, β, η) =α logE tr1 (ϑ; µ, Ω, α, β, η) = -d 2α + βt 2ηα 2 (t/αη) β-1 ,(68)D β logE tr1 (ϑ; µ, Ω, α, β, η) = D β logc(α, β) -1 2 (t/αη) β log(t/αη),

SPD matrix learning on synthetic dataset (Example 1). d 0 = 10, n train = 10 4 QRE 1.2×10 -9 5.8×10 -6 mMLP/ quad 2.3 × 10 -4 2.1 × 10 -2 MLP/ quad 2.5 × 10 -3 3.4 × 10 -1

Table content: <table><tr><td>Model/Loss</td><td>E quad</td><td>E QRE</td></tr><tr><td>mMLP/</td><td/><td/></tr></table>

Summary of the regression methods used in the experiments. = 200, d l = 20), Adam optimizer, activation functions as inTable 4 mMLP-mPE Section 4 3 layers (r l = 200, d l = 20), Adam optimizer, activation functions as in Table 4 MLP Goodfellow et al. (2016) 3 layers (200 units), ReLU activation, Adam optimizer, and default initialization of the scikit-learn 1 Gibbs iterations (discarded the first half), truncation of 5, and default initialization of the toolbox 2 . 1 scikit-learn: https://scikit-learn.org/ 2 NBCR: https://homes.cs.washington.edu/ ebfox/software/

Table content: <table><tr><td/><td>Method</td><td/></tr><tr><td>Acronym</td><td>Source</td><td/><td>Model specifications</td></tr><tr><td colspan="3">mMLP-mG 3 layers (r l GP Section 4 Rasmussen &amp; Williams (2006)</td><td>RBF kernel and default initialization of the sckit-learn 1</td></tr><tr><td>NBCR</td><td>Fox &amp; Dunson (2015)</td><td>10 5</td></tr></table>

Choice of activation functions for the mMLP-mPE and mMLP-mG models in Table3.

Table content: <table><tr><td>Model</td><td>µ</td><td>logη</td><td>Output layer α</td><td>β</td><td>Ω</td><td colspan="2">Hidden layers h H</td></tr><tr><td>mMLP-mG mMLP-mPE</td><td>linear() linear()</td><td>linear() linear()</td><td>-0.5+Sigmoid()</td><td>-0.5+Sigmoid()</td><td>* *</td><td>linear() linear()</td><td>* *</td></tr><tr><td colspan="2">* Trace-one Mercer Sigmoid:</td><td/><td/><td/><td/><td/><td/></tr></table>

Datasets and regression results.

Table content: <table><tr><td colspan="5">(a) Dataset specifications. Dataset</td><td/><td colspan="6">(b) Prediction errors in terms of the average RMSE across 10 runs. The bold-faced numbers indicate the statistical significance (p-value≤ 0.05).</td></tr><tr><td>Name oes10 edm atp1d atp7d scm1d scm20d</td><td>Appendix H.3.1 H.3.2 H.3.3 H.3.3 H.3.4 H.3.4</td><td>dx 298 2 411 411 280 61</td><td>dy 16 16 6 6 16 16</td><td>n train 50 50 100 33 1000 1000</td><td>n test 352 103 236 303 8802 8802</td><td>Dataset oes10 edm atp1d atp7d scm1d scm20d</td><td>mMLP-mPE 0.51 0.19 0.49 0.53 0.69 0.67</td><td>mMLP-mG 0.52 0.18 0.58 0.63 0.65 0.72</td><td>MLP 0.59 0.29 0.61 0.70 0.73 0.78</td><td>GP 0.97 0.32 0.89 0.97 0.82 0.87</td><td>NBCR 1.2 0.41 large large 1.5 1.7</td></tr></table>

Peer review (iclr2020):
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. While there is some interesting and potentially useful novelty in the approach, I have some concerns about the empirical evidence and modeling to truly determine the mMLP*s utility. First, in the synthetic data the mMLP with the l_QRE loss outperforms the l_quad loss with regards to both E_QRE and E_quad. Why does mMLP/l_QRE outperform on E_quad? As l_quad is actually designed to minimize this error, this is surprising and needs additional explanation. Additionally, there are two major changes between mMLP and the MLP. One is the network structure, and the second is the trace being normalized; which change really induces the improvements in the performance? If the l_QRE loss was used in the MLP, would you get similar improvements? The network structure as a whole needs greater validation. A major difference in (3) is that the weight matrices W_l are applied as a both a left and right multiplication. Given that the mathcal{H} operation symmetrizes and normalizes the matrix, a symmetric operation isn*t strictly necessary here. Using both multiplications leads to quadratic properties, which in my experience are less stable in the optimization. Can the authors validate this structure versus the simpler structure of simply using left multiplications? Or, in other words, is this weight multiplication structure helpful or do the benefits really come from the mathcal{H} operation. I think that the heteroscedastic regression experiments don*t evaluate on one of the key issues, which is uncertainty estimation. The trace-1 normalization is highly restrictive, so I imagine that that this method is getting the uncertainty incorrect. Also, heteroscedastic regression has a long history in neural networks, dating back to at least Nix and Wiegand in 1994. The manuscript needs to be updated to reflect the historical work and current literature on the topic. Please check Table 5(a), which states that you are only using a small number of training samples. Also, given the relatively small sample size of these datasets, please comment on the uncertainty of the results. How confident are you that the methods actually improve the prediction? How were the competing models tuned and optimized?