Matrix Multilayer Perceptron
Authors: Jalil Taghia
Maria Bånkestad
Fredrik Lindsten
Thomas Schön
Venue:  (2020)
Keywords: Multilayer Perceptron
symmetric positive definite
heteroscedastic regression
covariance estimation

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. 

Under review as a conference paper at ICLR 2020
MATRIX MULTILAYER PERCEPTRON
Anonymous authors
Paper under double-blind review
ABSTRACT
Models that output 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 (MIP). 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, introduce new variant of the MLP, referred to 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 scarce. We present an application of the model in heteroscedastic multivariate regression, including convincing performance on six real-world datasets.
INTRODUCTION
For certain applications, it is desirable to construct conditional covariance matrix 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-Fuclidean geometry (Pennec et al., 2005 ). The central problem that this work is concerned with 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 fully connected neural network (e.g..Goooll Goodfellow et al., (0166). 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-Fuclidean 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 such a way that it not only respects the SPD constraint, but also makes explicit use of it. We will achieve this by: 1) Explicitly taking the non-Euclidean geometry of the underlying SPD manifolds (c.g.. Pennee et al. 2005) into account by designing new loss function, and 2) by deriving 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 M The mMLP makes use of positive-defintt kernels to satisfy the SPD requirement across all layers. Hence, it provides a natural way of enabling deep SPD matrix learning.
We take step-by-step approach in the development of the model. We first develop 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 covariance regression setting (Section 3.2) application of the model discussed in the context of the heteroscedastic multivariate regression.
An implementation of this work will he made available via GitHuh.

---

Under review as a conference paper at ICLR 2020
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 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.. ongee 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 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 (lonescu et al. 2015 Huang & Gool 2017 Masci et al.) 20015] Huang et al., (2018). Our method falls in this category but differs in the problem Cormulation While previous methods address the problem where the input an SPD matrix and the output a vector, we consider the reverse problem where the input is a matrix with an arbitrary size and the output is an SPD matrix.
packpropan Our extension of the matrix backpropagation differs from the one introduced by Jonescu et al. (2015). In their work, the necessary partial derivatives are computed using 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 a-derivatives (Magnus) 2010) and its favorable generalization properties to derive a procedure which closely mimics the standard backpropagation.
MATRIX MULTILAYER PERCEPTRON
Matrix a-derivative. Throughout this work we adopt the narrow definition of the matrix derivatives known as the a-derivative (Magnus, 2010) in favor of the broad definition, the w-derivative. The reason for this is that the a-derivative has better generalization properties. This choice turned out to be crucial in the derivation of the mMI.P's backpropagation routine which involves derivatives of matrix functions w.r.t. the matrix of variables.
Preliminaries:
Definition: Let F be an m matrix function of an n x q matrix of variables X. The d-derivative of F(X) is defined as (Magnus 2010 Definition 2)
a ve)
DxF
a (vecX)TT
where F is mp nę 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 a-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 (ZI.... (Zi.....zz) .Z4) denote a matrix of variables z, Rd. The activation function XC(Z) defines a matrix function in the form of [K(Z)) = k(z,, 2,). Vi,j E (1......)) ..]....d....... where K is some differentiable activation function outputting scalar values. the following, we restrict ourselves to the kernel functions which form PD activation matrix functions. For numerical stability reasons, irrespective of the functional form K, we normalize the resulting matrix. This can be achieved by enforcing the trace-one constraint, HZ) = X(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 Mereer sigmoid kernel (Carrington al. (014) defined as
<(z,, 2,) tanh(oz; (3 @ tanh(az, + 1):
(2)
2

---

Under review as a conference paper at ICLR 2020
where and 3 denote the slope and the intercept, respectively. Furthermore, @ denotes the dot product. all experiments, we use default values of a = 1 and B = The The a-derivative of the Mereer sigmoid kernel is computed in Appendix E
Model construction. Let X E R Kp= indicate the input matrix and Y E Rdoxdo 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(Zn),
Zn = WoH, W + Bo.
H = H((z))
Z. WH WT + B,, AI L
(3)
H H((Z++1)
Z141 W 1. vee.v 1.+ 1 _ + +
The pair of W.ERdxxii v0 <I< L, and W RR are the weight matrices, B,E Rdexd, v0 50 < 1, are the bias matrices, Zr Rdd 40 <IS L+1, are the latent input matrices, and H2 E R4xdd VI < L 11, are latent output SPD matrices of trace one.
Design choice. In the construction of (3), we have ensured that H 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 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 T the normalized von Neumann divergence defined as:
AQRE + ( tr( tll@l $logg).
where log denotes the matrix logarithm (taking + as an example, it is computed using log VVdiagloog)) where V and are the matrix of eigenvectors and the vector of eigenval- ues 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 as
1 (AQRF YYI + AQRE( Y)).
EQRE Y, Y)
Taking a-derivative of QREE involves taking partial derivatives through the eigendecomposition. In the following, we derive method for analytically computing the derivative of ÉQRE-
The dederivative of the symmetrized von Neumann divergence. For the loss function defined in ( using the a-derivative's 's product rule and chain rule (refer to Appendix B). we obtain
! tret((( Yooog ! De Dettr(Yloggy 2
D l
(6)
where the above two terms are computed using
vec V - -YT) @ vec( a log
yY
vec VT - YT) @ vec( à
log Y)
De(tr((Y-- -Y)logY)
(vec(log Y) TTT
0Y.1
1xd8
vec T - -TT) o vec( a
logY)
&aando
1x
()(llYY))))) (vec(logY)T)
(8)
1xd3
Note that there multiple ways of symmetrizing the von Neumann divergence. Our choice in resembles the J-Jensen- --Jensen-Shhanon divergence for (Nielsen) 2010).
3

---

Under review as a conference paper at ICLR 2020
The remaining part in the computation of (6) is to evaluate
a logY, for all i, E aY.;
which involves taking derivatives through the eigendecomposition. In the following, we take a similar approach as Papadopoulo & Lourakis (2000 compute the necessary partial derivatives.
Let Y rdiag(A1..... Ado) Y T be the eigendecomposition. We can write
logY TAY
where A = diag(logAı. logAdo):
oY 0r Ar
(9)
+ T 2A r + TA
0r YY,
By multiplying (9) from left and right by r
and respectively, we obtain:
0r
0A
0r Y (
0A A (10)
r log T
A
+
,Y
0Y ij
0Y if
0Y,
0Y
where we have defined (2) T matrix, E, (T) + + ( = 0, which in turn follows from the fact that Y an orthonormal matrix,
Y and used the fact that B,y (Y) is an antisymmetric
oY.
Tr
Y Y OY | (2) ,, (Y) 0.
r Id
(11)
YY;
OY
Taking the antisymmetric property of E. (T) into account in I0]. we obtain
0Y log Ak= Ti TA
(12)
Biy (Tx) = Yix T Tr
(13)
2(logA: - logAs)
It is notable that by construction, do not have repeating eigenvalues, that is A F X Wk + there exists unique solution to (13 Once (T) is computed, it follows that
so
0T oi
or
YE (Y),
- (
(14)
oY
In summary, the necessary partial derivatives for computing (9) is given by 12 and computed for all j, we can evaluate (7) and ultimately evalaute (6).
14 Once
9 is
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 general version of the mMLP which produces both vector and SPD matrix as outputs. An important application of this model is in heteroscedastic multivariate regression which will discuss in Section 5.2
Model construction. As before, let X E xPxxP= denote the input matrix. The corresponding outputs in this case are: an SPD matrix of trace one E R x and y E R The mMLP of hidden layers is denoted by mMLP:X + { Y) and constructed as:
y b(zo),
2g C. YAyh, + bo:
Y H(Zo),
Zu = WoH, W + Bo.
h = b(z:),
2: = C,H,A + b
VI
(15)
H; H(Z)),
Z W,H W. +B,,
VI
h h(Z(+1),
Zr,+1 HI b
H = -H(Z++1)
Z = W (eW.veeX(WW++1pppp
BL+1.
where
h E R E IR dd
2r, b, E IR" Zr,
Z; B, E gdxd, Just as in the
C Ce v0 I+L+,, A RdR44 and W, E xdde41 VO <is L. standard MLP, h is an activation function of choice, e.g., the hyperbolic tangent function.
4

---

Under review as a conference paper at ICLR 2020
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
The remaining steps of feed-forward computation, backpropagation, and learning, are all described in Appendix D
Optimization.
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.
inputs x; E Rdd and a set of responses y; R [=1 be our training dataset consisting of a set of Consider the following multivariate regression problem: y(8) = fxx()) + (+)), where f is a nonlinear function, and e 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 Varly x.] for an unseen imput x, E Diest.
Model construction. Let Derain = (y(0), x (0))
We consider two likelihood models based Gaussian model and its generalization the multivariate power exponential model.
choices of the noise model, namely, the multivariate
Multivariate Gaussian model. We first consider the noise model follow zero-mean multivariate Gaussian (mG) distribution with a dense covariance matrix, that is e; ~ NO, Σ.). Let Σ; 77, 22, where 1/₃ tr(2 Ex) and 52; is trace-one matrix. The noise model can accordingly be reformulated as C; ~ N (0, S2; 17.) where Niri is referred to as the trace-one mG distribution Although these two formulations of the mG distribution are identical, we find it easier work with the latter. This is because the output layer of the mMLP model in (15) operates under the trace-one constraint. 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 easier to work with this reformulated mG distribution with trace-one covariance matrix which allows us to use the same choice of kernel function (2) across all layers.
Given e₁ (O, Dt my), the likelihood is defined
e == log Ni. (y: n4 S2, n),
mMLP
S2,
(16)
where Ely x] μl, and Varly x = n/22. The set includes all the neural network parameters.
Multivariate power exponential model. We next consider the noise model follow zero-mean multivariate power exponential (mPE) distribution (e.g., Gomez et al. (9988 which a generalized
variant of the mG distribution, that is (+) {0, E., oxe B,) where Σ, is symmetric real dispersion matrix and the pair of a; R+ and 3g E IR control the tail and the shape of the distribution. As special case, for a = 1 and [ = 1, the mPE includes the mG distributior
As in the case of the mG distribution, the dispersion matrix of trace-one. Let E, = n7,,, where 7/ = tr(E.). The noise model can accordingly be represented as ( ~ ELs (0, Ω2 04 Br m1) where Eut1 is referred to as the trace-one mPE distribution
find it easier to work with a reformulated variant where
Given ( Eur Ph,, C4i m)). the likelihood is defined as:
tem == log Eun1 (y: μr, 9 n). mMLP. x x {1,, a,
(17)
where Ely x] = and Varly x (av(3) where v(3) 2t/2T( d+1
The set @ includes all
the neural network parameters.
Optimization. Finally, the optimization involves learning the neural network parameters 0 by maximizing the likelihoods given by (16) and (17).
RRfeer to Appendix F.1 for the
functional form the trace-one mG distribution.
*Figure G.1 visualizes the probability density function of the distribution for selected values of a and 8.
'Refer Appendix G for the exact functional form of the trace-one mPE distribution and its hasic properties.
5

---

Under review as a conference paper at ICLR 2020
A. truth
mMLP, ÉQRE
mMLP, L Qquad
MLP, l Lquad
0.15
0:15
012
0.30
0.10
0.55
0.0
585
0.0
0.00
5.00
-5.3
000
-0.10
-2.0
0.2
0.12
0.30
9.05
5.02
0.05
0.00
-0.0
00.00
-8.35
-9.05
-0.10
Figure 1: SPD matrix learning (refer to Example in Section (.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 fquad
EXPERIMENTS
Experiments are divided into two parts. The first part is on the empirical validation of the mMLP model in supervised task of learning SPD matrices using synthetic data. The second part discusses our results in heteroscedastic multivariate regression on real-world datasets.
5.1 SPD MATRIX LEARNING
5.1.1 THE CHOICE OF Loss FUNCTION
Consider the problem of learning SPD matrices on synthetic data using the mMLP model of (3). The objectives are 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 normalized von Neumann divergence fORE Y Y) given by (5). The LQrE is related to the Riemannian geometry. The second candidate is the quadratic loss, lquad Y,, Y) tr((( Y) which related to the Euclidean geometry.
Example 1. Consider the set Dtrain= (x, Y, ) M of inputs x; E To:20 and corresponding SPD matrix outputs Y, d which are in this case dense covariance matrices (refer to Appendix H.1 for details on the data generation). The goal to estimate the covariance matrices Y associated to the input vectors from the unseen test set Dtest {x J The training size is varied between Terain = (10?, 10*} samples. The analysis is carried out for do = {10, 20}. Two examples of the test target outputs Y, for dg 20 and Ptrain 102 are visualized in Figure I
The mMLP models 3 are trained using layers (20 units per layer, d, = 20) under our two choices of loss functions: εQRE, and fquad All models share the initialization and the only difference here is the loss function. Refer to Appendix [or additional details on the mMLP initialization the test set, Nlasst 10 in terms of the losses the error measures, shown as FORE and Equad Table [ssmmarizes the results of the evaluation (also refer Figure 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 {oRE that takes into account the geometry of the SPD matrices clearly outperforms the one based on the Euclidean geometry, qqaud- The advantage is more pronounced for small set of training data.
The performance evaluated
5.1.2 THE MMLP VS. THE MLP BASED ON CHOLESKY DECOMPOSITION
As discussed in Section ! one can take heuristic approach and tailor the output layer of vanilla MLP so that the resulting output matrix satisfies the SPD requirement.
For the sake of comparison, we solve the same task in Example I using the standard MLP with layers (200 units per layer) and with the quadratic loss. To meet the SPD requirement, we use the
6 'Her and in general throughout this section, special care has been made to minimize the effect of overfitting through trying various initializations and using early stopping.

---

Under review as a conference paper at ICLR 2020
Table 1: SPD matrix learning on synthetic dataset (Example 1).
dg
Pirain
10
de 20, Plrsinn
10
Parain
102
Perain 10
ModeVT
Rquad
FQRE
Equad
EquE
Equad
EqnB
Equad
EquE
MMLP/EQREE
1.2x10
6.8x10
3.3x10-a
1.7x10-5
6.3x10-6 11x10-4
4.3x10-5
222x10
mMLP/
2.3 x 10
2.1 10
2.1 10-3
3.9 10
3.7 10
6.6 10
6.0 10
1.2
MLP //quudd
2.5 10-3
3.4x
13x
811
6.5
4.6
6.4
Table 2: Shallow vs deep architecture, Example 2. Error is measured in terms of EQRF.
do
Model
2x 6 L-
7
L-4
1=2 L=4 L-6
Deep architecture
B
3 3x110 4 46x10 4x10 8x10-2 6x10 7XX1-
Shallow architecture
18
4x10-
6X10-3
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 unique diagonal matrix with positive diagonal entires. Taale|]summmnt quite poor for small set of training data. As the size of training data grows, the performance improves as expected (refer to Figure D for the visual illustration).
the results of the evaluation. Overall the performance is
5.1.3 SHALLOW DEEP SPD MATRIX LEARNING
The design of the mMLP model 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(Zo),
Zo Woh (Wol) + Bo:
h: b(2)),
2j Whh by, 1 1515L,
(18)
h = ((2++1),
L W W + bL+1-
This amounts to 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 captures complex dependencies which are important for the SPD matrix learning. Below, present 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 for the case of Merain 10 and output dimensions do {10, 20). We directly compare the performance of 33) against (18) under different number of hidden layers = (2, 4, 6). For the mMLP model, the number of hidden units at each layer is set to 20, and for the MLP model is set to 200 units. The shallow design (18) uses the hyperbolic tangent as the activation function (-). The same choice of the activation matrix function H(-),, given by is used for both models. We use ÉQRE as the choice of the loss function for both models (refer to Appendix H.2 additional details on the initialization). The performance is evaluated terms of EQRE.
Table summarizes the results of the evaluation. Although the shallow design 18 performs relatively well, underperforms comparison (3). Given the limited number 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 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 model can effectively capture the heteroscedasticitty then it will provide better uncertainty quantification, and the estimation the conditional mean response should also improve, in comparison to the model that build based the homoscedasticity assumption.
7

---

Under review as a conference paper at ICLR 2020
Table 3: Summary of the regression methods used in the experiments.
Method
Acronyus
Source
Model specifications
P-mG
Sectioe
layers (fr 200, de 200, Adam optimizer, activation functions
Sectwod -
layers (F;
200, de 20), Adamn optimizer, activation functions labley
MI.P
GoodfflTww al. 2016
layors
Rel.U octivation,
optimizer, and default initializztion the
GP
Ravmmosen William
2006
RBF kerol and default intisaizzation of the sckit-leart"
NBCR
Fox & Dunson
2015)
10" Gibhs Herations (discanted the fist hair), truncction of 6. and defonh initaizatiion of the soolbor
seikät-lramm: tripss/sskkiikeenn..c ong
2
NBCR: https//liouess washingto edut ohfouhotwarre
Table 4: Choice of activation functions for the mMLP-mPE and mML.P-GG models in Table 3
Outpul bayer
Hidden layers
Model MM.PParii
loga
f2
(
lnear ()
linear()
IMLP--PE
linear()
linear()
0.64Sigmoid()
0.5 0.5- Sigmoid()
linear()
Trace-une Mercer Sigmoid: (2))
K (2)
where (K(Z)),,) = n( n(),,,)), V5,J. J. The kemel function x(-, given by
Table 5: Datasets and regression results.
(a) Dataset specifications.
(b) Prediction errors terms of the average RMSE
across 10 runs. The bold-faced numbers indicate the statistical significance (p-value 0.05).
Dataset
Name
Appendix
dx
dy
Plerain
Vitest
Dinset
mMI P-mPF
n-M.P-mfi
MIP
GP
NBCR
oes
1.3.1
298
16
50
352
0.51
edm
-30
103
edm
0.19
0.18
0.29
0.32
0.41
splid
apidd
0.49
0.58
0.61
0.89
large
up7d
411
33
303
0.53
0.70
0:97
semld
280
16
1000
8802
somld
smm20d
1000
8802
scm20d
0.67
0.72
0.78
0.87
For this purpose, we compare our mMLP- regression models against another heteroscedasti- based covarccoe regression model by Fox Dunson (2015) and two homoscedastic- mean-regression models, namely, the MLP regressor and the Gaussian process (GP). The models used this experiment are summarized in Table
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 sets, shown in Table 5b The result suggests that the mMLP homoscedastic-basee regression models are capable of capturing dependencies between the output measurements which contributes to the better estimation of the mean predictions.
error (RMSE) on test
LIMITATIONS AND FUTURE WORK
The main limitation the mMLP has to do with scalability to higher dimensions. The complexity associated with computing the a-derivative of the von Neumann loss function (5) at the output layer is (dg). Taking the symmetric nature of the SPD matrices into account, the computational complexity the hidden layer i reduces to (df)).
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 (df x df). However, these matrices are in fact available in the form of sparse block matrices, which means that it possible 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 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.
DISCUSSION
We introduced new method to lear SPD matrices, referred to the matrix multilayer perceptron (mMLP). The mMLP takes the non-Euclidean geometry of the underlying SPD manifolds into account by making
of the von Neumann divergence the choice of the SPD manifold metric. One key aspect of the mML.P 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 heteroccedastic multivariate regression. We showed the effectiveness of the proposed model on multiple real-world datasets.
8

---

Under review as a conference paper at ICLR 2020
REFERENCES
Vincent Arsigny, Pierre Fillard, Xavier Pennec, and Nicholas Ayache. Geometric means in novel vector space structure on symmetric positive-definite matrices. SIAM Journal of Matrix Analysis Applications, (9(1)328-344) 2006.
André M. Carrington, Paul W. Fieguth, and Helen H. Chen. new Mercer sigmoid kernel for clinical data classification. In International Conference of the IEEE Engineering in Medicine and Biology Society, 2014.
Mohammed E. Fathy, Azadeh Alavi, and Rama Chellappa. Discriminative log-Euclidean feature learning for sparse representation-basee recognition of faces from videos. In I/CAI, 2016.
Emily B. Fox and David Dunson. Bayesian nonparametric covariance regression. Journal of Machine Learning Research, 16(12):2501-23) 2015.
Gómez, Miguel Gomez- Gomez-Villegas, and Marin. A multivariate generalization of the power expo- nential family of distributions. Communications in Statistics-theory and Methods, 27(3):589-600 1998.
Ian Goodfellow, Yoshua Bengio, and Aaron Courville. Deep Learning. MIT Press, 2016
Mehrtash T. Harandi, Conrad Sanderson, Richard Hartley, and Brian C. Lovell. Sparse coding and dictionary learning for symmetric positive definite matrices: A kernel approach. In ECCV, 2012.
Mehrtash T. Harandi, Mathieu Salzmann, and Richard Hartley. From manifold to manifold: Geometry- aware dimensionality reduction for SPD matrices. ECCV, 2014.
Peter D. Hoff and Xiaoyue Niu. covariance regression model. Statistica Sinica, 22(2))729-753, 2012.
Zhiwu Huang and Luc Van Gool. A Riemannian network for SPD matrix learning. AAAI, 2017.
Zhiwu Huang, Jiqing Wu, and Lue Van Gool. Building deep networks on Grassmann manifolds, In AAAI 2018.
Catalin Ionescu, Orestis Vantzos, and Cristian Sminchisescu. Matrix backpropagation for deep networks with structured layers. In ICCV, 2015.
A. Karalic and I. Bratko. First order regression. Machine Learning, 26(2-3):147-176, 1997.
Diederik P. Kingma and Max Welling. Auto-encoding variational Bayes. In ICLR, 2014.
Jan R. Magnus. On the concept of matrix derivative. Journal of Multivariate Analysis, 101(9): 2200-2206, 2010.
Jonathan Masci, Davide Boscaini, Michael M. Bronstein, and Pierre Vandergheynst. Geodesic convolutional neural networks on Riemannian manifolds. In ICCVW, 2015.
Ha Quang Minh, Marco San-Biagio and Vittorio Murino. Log-Hiibert-Schmid metric between positive definite operators on Hilbert spaces. In NeurIPS, 2014.
Frank Nielsen. A family of statistical symmetric divergences based on Jensen's 's inequality. arXiv, 2010.
Frank Nielsen and Richard Nock, On the centroids of symmetrized Bregman divergences. CORR, 2007.
Michael A. Nielsen and Isaac L. Chuang. Quantum Compulation and Quantum Information. Cam- bridge University Press, 2000.
Peter Meer Oncel Tuzel, Fatih Porikli. Pedestrian detection via classification on Riemannian man- ifolds. IEEE Transactions on Pattern Analysis and Machine Intelligence, 30(10):1713-1727, 2008.

---

Under review as a conference paper at ICLR 2020
Théodore Papadopoulo and Manolis I. A. Lourakis. Estimating the Jacobian of the singular value decomposition: Theory and applications. In ECCV, 2000.
Xavier Pennec, Pierre Fillard, and Nicholas Ayache. A Riemannian framework for tensor computing. International Journal Computer Vision, 66(1):41-66 2005.
Mohsen Pourahmadi. Joint mean-covariance models with applications longitudinal data: Uncon- strained parameterisation. Biometrika, 86(3):677-690, 1999.
C. E. Rasmussen and C. K. I. Williams. Gaussian processes for machine learning. New York, NY, USA, 2006.
David E. Rumelhart, Geoffrey E. Hinton, and Ronald J. Williams. Learning representations by back-propagating errors. Nature, 323(10):533-53, 1986.
Eleftherios cerromitros-iioufis Grigorios Tsoumakas, William Groves, and Ioannis Vlahavas. Multi-target regression via input space expansion: treating targets as inputs. Machine Learning. 104(1):55-995 2016.
Koji Tsuda, Gunnar Rätsch, and Manfred K. Warmuth. Matrix exponentiated gradient updates for on-line learning and Bregman projection. Journal of Machine Learning Research, 6(1):995-1018, 2005.
10

---

Under review as a conference paper at ICLR 2020
A MATRIX NOTATIONS
We use T for the transpose operator, tr(-) ) for the trace operator, and det(.) for the matrix determinant. The symmetric part of square matrix B is denoted by sym (B) (B B //2. The Kronecker prod- uct is denoted by &, the Hadamard product by 9, and the dot product by O. We use the vec-operator column-by-column stacking of a matrix A, shown as vec = vec A). Let be an m n matrix, the operator P (m..) will then rearrange vecA to its matrix form: A = (m.n) (vecA). For the m X n dimensional matrix A, the commutation matrix is shown as K (m,m) which is the mn x mn matrix that transforms vecA into vecA T as: K (m..) = vecA An m X m identity matrix is shown as m. If T(X:: Raxd - R is real-valued function on matrices, then VxT(X) denotes 0T
the gradient with respect to the matrix X, Vx7(X) = 884 i,j=lad' The matrix logarithm and the matrix exponential are written logA and A, respectively. The matrix exponential in the case of symmetric matrices can be expressed using the eigenvalue decomposition as epp = (errp "(TT V where is an orthonormal matrix of cigenvectors and A 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 symmetric and strictly positive definite (PD) it computed using logA V (logA)V where (IogA)i,: logAi,i
B THE a-DERIVATIVE: DEFINITION AND PROPERTIES
Definition Let F be matrix function of q matrix of variables X. The a-derivative of F(X) is defined as (Magnus 2010, Definition 2)
Dx = vecF
(19)
δ (vecX)
where Fx is an mp ng matrix which contains all the partial derivatives such that each row contains the partial derivatives of one function with respect to variables, and each column contains the partial derivatives of all functions with respect to one variable.
Let F (m x p) and G (p x r) be functions of X (n q). Then the product rule for the a-derivative is given by [Magnus] 2010)
Product rule
DxGF = (G @ Lm)DxFFF + (1, & G.
(20)
Chain rule Let F (m x p) be differentiable at X (n x ②) and G ( x r) be differentiable at Y F(X), then the composite function H G(F(X)) is differentiable at X, and
Dx GDGFF
(21)
which expresses the chain rule for the a-derivative MMagnus 2010
C THE BASIC CASE OF THE MMLP
C.I FEEDFORWARD STEP
At the feed-forward computation, we compute and store the latent outputs Y, H. for all 1,....,)
I E (T+1,, 1,..,, 1) using the current setting of the parameters, which are W,, and B, computed from the learning step, Appendix C.3
C.2 BACKPROPAGATION STEP
We first summarize the necessary a-derivatives for the backpropagation, and then write down the backpropagation procedure accordingly.
Derivatives required for backpropagation The derivative of the activation matrix function de- pends the specific choice of kernel function, and in general it is computed readily from the

---

Under review as a conference paper at ICLR 2020
definition of a-derivative,
(Zi)((z))
Dz, HZ,)) = a (vecZi) )
LE ++1...
(22)
For our specific choice of activation function, the Mercer Sigmoid kernel Appendix H
2 it is computed in
Via repeated use of the product rule of a-derivatives (20) we obtain
Dw, Zi (Wr @ (H @ La) (Ld, 8 (W)H++1)K (di,dit))) E (O....,,4 4}, (23)
Dw.. Z, = (W 5+ 1 pipe @ d ((vecX) & Idrt1)
(I @ (W WI++ vec d2 @ ) (dz 1.112))
(24)
DH. 2 = (W, @ (1 @ Wi), E LJ,
(25)
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 Hg Y. Backpropagation to the hidden layer is computed recursively using the derivatives computed at the previous layer according to
&z, (z, Hi,
v/ {0, (O....,, L 1},
(26)
Dw. l Dz, &Dw Zi,
WE {0, ......... + 1}, v (0..... EE (0,...,L+11,
(27)
DH { Dz, (DH.,, Zt.
(28)
DB, l Dz, L
(29)
C.3 LEARNING STEP
Learning involves updating the weights W, and the biases B, using derivatives computed during backpropagation step. These are updated using derivatives Dw, and DB, for given learning rate
as
W. + W, r/P (di de+ (Dw,e), E)
W (0, (O...., L 1},
(30)
B, + B (d., d:) (DB, E),
W E (0, 1},
(31)
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, ÿ and H,, Y using the current settings of the parameters, for all / (0, L+
D.2 BACKPROPAGATION STEP
Most of the necessary derivatives are identical the ones in Appendix C.2 However, there 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 computed using the definition of the a-derivative,
0 b
Dz. bi(z))
LE (0, ...., L 1}.
(32)
(zi) T'
12

---

Under review as a conference paper at ICLR 2020
The other required derivatives computed
Z C,H,
LE [O,...,)),
(33)
DA 2L+1 = (CL+1 (CL+)H @ Idmm)):
(34)
D (H: T & I
LE (0,..... D),
(35)
Dc A4+1 = (HL+1A A & Ium
(36)
Di Zi C,
Le {O..... L}
(37)
DH, = ((Arh++1) & (Ia & Ci),
Le {O....,, L}
(38)
D 21+1 = ((A4+1let)) & (rdr) & C1+1).
(39)
Backpropagation For simplicity of expressions, let ho ÿ and Ho Y. The derivatives recursively computed as
Dng f = Dye
(40)
DHo Dy = D20 (Dz +
Del
(41)
2z l Dr. (D, hi,
W {0, ........................ 1},
(42)
D = Dz, (D
E (0........................ L
(43)
Dz, / = DH, (Dz, HI,
W
+ 1},
(44)
DH, { 22111 (DH, Ze+1 + (D,,, Zi
A1 E (0. (0.....4) L))
(45)
DA, DR. 2,
W E (0,
1+1},
(46)
Dc, l = Dz, &D Z1,
VI E (
(47)
Dw, l Dz, DDw Zi.
A E (0,
L+1),
(48)
Db, Dz. e
M E (0.....
(49)
DB, Dz, L,
W E (0, 11.
(50)
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,, and B, are the same as the ones given in Appendix
C.3 The remaining parameters are learned in a similar fashion,
A A r/P (d. r (D)
VI E (0, ... L + 1},
(51)
C C T/P (mid) (Dc,
V/ E (0, .., L 1},
(52)
b + b 7/P (m
VI E (O...., L 1}.
(53)
E THE a-DERIVATIVE OF THE MERCER SIGMOID KERNEL
The a-derivative of the Mercer sigmoid kernel can be computed as
Kmn
DzH
Oz, trk
Vi,m,n .l.....dij.
(54)
1x
where z; indicates the 7th column of Z, trK = triC(Z), Kmn = <K(mm 2n) as defined in and
i Eynn
(a(1T - f(z,) (2ii))
(trCC
(Zm 2n)
22imm 22 (55)
8z, trk
(trZ)²
Ozy
where f(z;) tanh(az, B)
13

---

Under review as a conference paper at ICLR 2020
10.1, 0.010
[0.01,
19.0%,
(0.01.2.))
(0.01 10.00
7.5
430
63
43
616
45
16
33
0.00)
11.0, 0.5)
(1.4, 1.20
(1002.00)
(1.0.10.0)
0.25
0.30
4.2
0.893
0.06
4.5
0.45
0.16
23
044
0.093
0.12
62
0.00
0.15
0.993
6/0
1e:
0.511
5.07
0.899
4.44
2.08
4.40
5.06
0.39
0.39
0.37
6:20
0.44
0.45
0.355
0.238
0.06
0215
Figure G.1: The probability density function of trace-one mPE distribution, (1, M/ 1J, 0,
for fixed values of p. 1) tr(MM and varying values of scale and shape parameters (a,, 8). When and B = 1, the density corresponds to the multivariate Gaussian distribution. N(u, M).
F TRACE-ONE MULTIVARIATE GAUSSIAN DISTRIBUTION
F.1 PROBABILITY DENSITY FUNCTION
For d-dimensional random variable Rd, we define the trace-one Gaussian distribution according
to
01 del (2mm)22)
(0- ((n2)))
(56)
where μe E IRd is the mean, 7) E R+ is the scale parameter, and f2 E: Fdxd x is the trace-one covariance matrix, tr(22 1.
F.2 THE a- A--EEIVATIVES OF THE L.OGPDD
The a-derivatives the trace-one Gaussian distribution' 's log--dd with respect to its parameters are summarized
Dn (1: - (vec((2 1))
(Dnt)
(57)
Dpl = -ec - m)(0 - 40 12
(58)
Dn Ω, 11) (0 (1((112)-11
(59)
1
Dong = 2 + (0 μ1) (n))) (1 - pd.
(60)
2
G TRACE-ONE MULTIVARIATE POWER EXPONENTIAL (MPE) DISTRIBUTION
G. PROBABILITY DENSITY FUNCTION
For a random variable 0 E TRd the trace-one mPE distribution can be expressed as
c(a,
1 1/t(0::..0 2252)
12,
exp
()e(())
2
(61)
(T
da, 3)
100;u, = (0 41) T S (0 4), tr((2)
1.
n T
Here, 4 is the mean vector, and Ω is a x d symmetric real dispersion matrix where tr( 1. The parameter 7 has the same role in the trace-one Gaussian distribution. The pair of R
14

---

Under review as a conference paper at ICLR 2020
and [ E RR* control the tail and the shape of the distribution. As special case, the mPE includes the Gaussian distribution: For a = and 3 = 1, the trace-one mPE distribution corresponds to the trace-one multivariate Gaussian distribution. Figure G visualizes the probability density function of the distribution for selected values of a and B, for the case of d 2.
Remark: Very large and very small values of and 3 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 bound them conservatively as: 0.5 S a, S 1.5, using the sigmoid function.
G.2 MOMENTS
Let 0 E Rd and &- (H, a, B.n). puted from:
The mPE's mean vector and covariance matrix are com-
2 412)
0 μ
V|ø
(3)(,
((8)=
(62)
dl 23
where I'()) denotes the gamma function.
G.3 THE a-DERIVATIVES OF THE LOG-PDF
It is straightforward to take derivatives of the mPE's log-pdf using the favorable generalization properties of the a-derivative's chain and product rules. These are summarized as:
D (0; u, S2 a, 3, 11) - (vec (ve
B
(t/an)*-1Dnt
(63)
2001
Dat (vee(S1 10 u)(0 u) =
(64)
B (1/an))-1
D S2, = 20m
(65)
Du 11) TT)
(66)
d
Bt
D DoggllogEu 14, =
+
(1/om) 8-1
(67)
2
207
(0: 3,11 d
3t
(1/00)
(68)
2a
27a2
DglogEir (0: 44 S, B,n) = Dzlogc(a, B)
(t/ (t/om) log(t/om))
(69)
2
1
Dølogo(a, B) 3 2,32 (v(d)23) + log2).
(70)
H ADDITIONAL DETAILS ON THE EXPERIMENTS
H.1 EXAMPLE
Data generation Let A ddo be matrix where each of its elements is generated from standard normal distribution. The matrix A is kept fixed. The 2th class covariance Y, computed according to the following procedure:
1. Draw 10 samples from a known Gaussian distribution ((,,, E,) with unique mean f4 TR and unique dense covariance matrix Σ, E Ro
2. Let t, be a random sample from this Gaussian. For this sample, compute y, = Aty. For all
10 samples, collect the results into = {o,
3. Compute the sample covariance of y and normalize the resulting covariance matrix to trace one, that is Y, (ow(y/tt(oon)
15

---

Under review as a conference paper at ICLR 2020
Initialization All models use the same batch size (equal 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).
H.2 EXAMPLE 2
Data generation
See the data generation procedure in Example 1.
Initialization Both models 18 and use the same batch size (equal to 5), the same choice of loss function ( 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 the activation function across the hidden layers, while (3) makes use of the same choice of the activation matrix function in its output layer.
H.3 MULTI-OUTPUT REGRESSION DATASETS
H.3.1
oes10
The dataset ocs10 was obtained from Spyromitros-Xione Xioutiss al. (2016). The Occupational Em- ployment 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-Xion Xioufis et al. 2016), "each row provides the estimated number of full-time equivalent employees across many employment types for 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 Table Sa).
H.3.2 edm
The dataset edm obtained from Karalic & Bratko (997).. 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 Sa
H.3.3 atpld AND atp7d
The datasets atpld and atp7d were obtained from (Spyromitros-Xinne Xioufis et al., 2016). The Air- line Ticket Price (ATP) dataset includes the prediction of airline ticket prices. As described in "Spyromitros-Xiouff et al.] 2016), the target variables are either the next day price, atpld, or minimum price observed over the next seven days atp7d for target flight preferences listed in Spyrroitross Xiouñi 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-weck of the observation date. Data samples were randomly divided into training and test sets (refer to Table 5a).
H.3.4 semld AND scm20d
The datasets scmld 1d and scm20d were obtained from (Spyromiros-Xioutis 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 specific tournament day. The datasets contain 16 regression targets, where each target corresponds to the next day mean price scmld or mean price for 20 days in the future scm20d for each product Sppromitros-Xioutt et al. 2016 Table 5). Data samples were randomly divided into training and test sets (refer to Table 5a).
16


Reviews:
[{"review_id": "ByxIxKl6YS", "reviewer": "AnonReviewer2", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper generalized neural networks into case where a semidefinite positive matrix is learned at the output. The paper presents theoretical derivations that look sound, and validating experiments on synthetic and real data. I must say my expertise does not really correspond to what is done in this paper, but I do not see any obvious flaws and the results look solid. I appreciated the discussion of limitations in section 6. I vote for acceptance with the weakest possible confidence level since it is likely I missed many important points. ", "extra_scores": {"experience_assessment": "I do not know much about this area.", "review_assessment:_thoroughness_in_paper_reading": "I made a quick assessment of this paper.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I did not assess the derivations or theory."}, "extra_text": {}}, {"review_id": "BkgR_zORKS", "reviewer": "AnonReviewer1", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "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. \n\nPros: \n1. The use of von Neumann divergence as a loss for this task is perhaps novel.\n2. The use of \\alpha-derivatives, while computationally demanding, is perhaps novel in this context as well. \n\nCons:\n1. I do not think the problem setting or the proposed framework is entirely new or is the best choice of its ingredients. Specifically, the idea of using second-order neural networks have been attempted in several prior papers, including the ones the paper cite (such as Ionescu et al. ,2015). Several other papers in this regard are listed below. \n[a] Second-order Temporal pooling for action recognition, Cherian and Gould, IJCV, 2018\n[b] Deep manifold-to-manifold transforming network, Zhang et al, ICIP, 2018\n[c] Statistically motivated second-order pooling, Yu and Salzmann, ECCV, 2018\n\nIn comparison to these methods, it is not clear how the proposed setup is novel, or in what way is method better. There are no comparisons to these methods, and thus it is difficult to judge the benefits even empirically. \n\n2. There are also models that predict the mean and covariance matrices directly from the model, such as the works below. The paper should also include and perhaps compare to their datasets.\n[d] Deep Inference for Covariance Estimation: Learning Gaussian Noise Models for State Estimation, Liu et al, ICRA, 2018\n\n3. I do not think the use of von Neumann divergence as a loss is the best choice one could have, esp. for a deep neural network learning setting. This divergence includes the matrix logarithm, which is perhaps computationally expensive. This is a problem when using other popular loss/similarity functions on SPD matrices (such as the log-Euclidean metric). Perhaps a better option is to use the Jensen-Bregman log-det divergence, as suggested in [a] above; this divergence is symmetric and also has computationally efficient gradients. It is unclear why the paper decided to use von Neumann. \n\n4. The experiments are not compelling, there are no comparisons to alternative models and the datasets used are small scale. Thus, it is unclear if the design choices in the paper have any strong bearing in the empirical performances.\n\nOverall, the paper makes an attempt at designing neural networks for learning SPD matrices. While, there are some components in the model that are perhaps new, the paper lacks any justifications for their choices, and as such these choices seem inferior to alternatives that have been proposed earlier. Also, the experimental results are not convincing against prior works. \n\n", "extra_scores": {"experience_assessment": "I have published in this field for several years.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I assessed the sensibility of the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I carefully checked the derivations and theory."}, "extra_text": {}}, {"review_id": "BJg1Bzy0Yr", "reviewer": "AnonReviewer3", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "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.\n\nWhile 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.\n\nFirst, 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?\n\nThe 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.\n\nI 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.\n\nPlease 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?", "extra_scores": {"experience_assessment": "I have published in this field for several years.", "review_assessment:_thoroughness_in_paper_reading": "I read the paper thoroughly.", "review_assessment:_checking_correctness_of_experiments": "I carefully checked the experiments.", "review_assessment:_checking_correctness_of_derivations_and_theory": "I assessed the sensibility of the derivations and theory."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper introduces a novel architecture and loss for estimating PSD matrices using neural networks.  There is some theoretical justification for the architecture, and a small-scale but encouraging experiment.

Overall, I think there is a sensible contribution here, but there are so many architectural and computational choices presented together at once that it's hard to tell what the important parts are.

The main problems with this paper are:
1) The scalability of the approach O(N^3)
2) The derivation of the architecture and gradient computations wasn't clear about what choices were available and why.  Several alternative choices were mentioned but not evaluated.  I think the authors also need to improve their understanding of automatic differentiation.  Backprop through eigendecomposition is already available in most autodiff packages.  It was claimed that a certain kind of matrix derivative provided better generalization, which seems like a strong claim to make in general.
3) The experimental setup seemed contrived, except for the heteroskedastic regression experiments, which lacked competitive baselines.  Why were the GP and MLPs homoskedastic?

As a matter of personal preference, I found that having 4 different "H"s differing only in font and capitalization for the network architecture was hard to keep track of.

I agree that R1 had some unjustified comments and R2's review was contentless.  I apologize for these inadequate reviews. 
Decision: Reject