where $\mathbf{1}_{q}$ is $q$ dimensional vector of 1’s. In this particular
case, the latent variable $z$ can be alternatively represented as
$z=B^{\mathrm{\scriptscriptstyle{T}}}x+\rho\mathbf{1}_{q}u+e$ where
$u\sim{\mathrm{N}}(0,1)$ and $e\sim{\mathrm{N}}(0,1-\rho)$ independent of $u$.
Hence conditional on $u$, univariate probit models suited to INLA can be
fitted with the augmented covariates $\tilde{x}=[x,u]$ to each column of $y$.
Furthermore, if $\tilde{\beta}_{j}$ represents the $j$th column of the
regression coefficient matrix $\tilde{B}$ for this conditional probit model,
then an estimate of $\beta_{j}$ can be formed by setting
$\hat{\beta}_{j}=(\sqrt{1-\rho^{*}})\tilde{\beta}_{j}$. We implement this
along with a random walk MH algorithm for updating $u$ where the standard
deviation of the proposal distribution is updated following Haario et al.
(2001) so that the acceptance rate is optimized to $\sim 0.23$.

We conducted two separate experiments with the true correlation
$\rho^{*}=0.3,0.5$ and with each $B^{*}$ generated according to the “Dense”
case scenario described in Section 4 of the main document. For each case we
fix the sample size $n=200$, $p=3$ and consider $q=3,5,7,9$; we did not
consider higher dimensional responses as the INLA-MCMC method did not scale
well computationally with growing $q$. We report the errors in estimating the
regression coefficients, namely $\|B^{*}-\hat{B}\|_{F}^{2}/q^{2}$, averaged
over 30 independent replications where the prior distribution of outcome-
specific regression coefficients are ${\mathrm{N}}(0,5^{2})$ and an LKJ(1)
prior on the correlation matrix so that the marginal priors on the
correlations are $\text{Beta}(q/2,q/2)$. When executing the rescaling of the
coefficients we set $\rho=\rho^{*}$. We summarize the results in Table S.1.
The results show that bigMVP performs better than INLA-MCMC in all the cases
considered here although it should be kept in mind that the rescaling of the
estimates from INLA-MCMC were computed assuming the truth is known.

| $q=3$ | $q=5$ | $q=7$ | $q=9$
---|---|---|---|---
| bigMVP | INLA-MCMC | bigMVP | INLA-MCMC | bigMVP | INLA-MCMC | bigMVP | INLA-MCMC
$\rho^{*}=0.3$ | 0.04 | 0.07 | 0.06 | 0.09 | 0.07 | 0.11 | 0.04 | 0.08
$\rho^{*}=0.5$ | 0.05 | 0.08 | 0.05 | 0.09 | 0.06 | 0.10 | 0.04 | 0.09

Table S.1: Comparison between bigMVP and INLA-MCMC in estimating the
regression coefficients when the correlation structure between the binary
outcomes are of the form
$\Sigma=(1-\rho^{*})\mathrm{I}_{q}+\rho^{*}\mathbf{1}_{q}\mathbf{1}_{q}^{\mathrm{\scriptscriptstyle{T}}}$.

##  S.6 Marginal pairwise prediction

Although our focus is on inference on the model parameters, it can nonetheless
be useful to obtain predictive distributions as key components of model
assessments and comparisons. As obtaining a joint predictive rule for all $q$
outcomes is necessarily computationally demanding when $q$ is large, we focus
on prediction of an arbitrary pair of outcomes in a new sample $t$ having
covariate value $x_{t}$.

Fix a pair $(j,k)$. Then from our proposed two-stage approach, we have access
to the approximate marginal posterior distributions $\Pi_{j}^{*}(\beta_{j}\mid
y,X),\Pi_{k}^{*}(\beta_{k}\mid y,X)$ and $\Pi_{jk}^{*}(\sigma_{jk}\mid y,X)$.
Let
$\theta_{jk}=\\{\beta_{j}^{\mathrm{\scriptscriptstyle{T}}},\beta_{k}^{\mathrm{\scriptscriptstyle{T}}},vec(\Sigma_{jk})^{\mathrm{\scriptscriptstyle{T}}}\\}$.
Given a new test point $x_{t}$, the marginal predictive distribution for this
pair is

$p_{jk}(y_{t}^{(j)},y_{t}^{(k)}\mid y,X)=\int
p_{jk}(y_{t}^{(j)},y_{t}^{(k)}\mid\theta_{jk},y,X)\Pi(\theta_{jk}\mid
y,X)d\theta_{jk},$ (S.12)

where we approximate the joint posterior $\Pi(\theta_{jk}\mid y,X)$ by a
product of the marginals, $\Pi(\theta_{jk}\mid
y,X)\approx\Pi_{j}^{*}(\beta_{j}\mid y,X)\Pi_{k}^{*}(\beta_{k}\mid
y,X)\Pi_{jk}^{*}(\sigma_{jk}\mid y,X)$. The first component inside the
integral in (S.12) can be written as
$p_{jk}(y_{t}^{(j)},y_{t}^{(k)}\mid\theta_{jk},y,X)=\int
p_{jk}(y_{t}^{(j)},y_{t}^{(k)}\mid
z_{t}^{(j)},z_{t}^{(k)},\theta_{jk},y,X)dz_{t}^{(j)}dz_{t}^{(k)}$, where
$(z_{t}^{(j)},z_{t}^{(k)})$ are the corresponding latent variables for the
test point $x_{t}$. Hence, we focus on the predictive distribution of the
latent variables $(z_{t}^{(j)},z_{t}^{(k)})$. Under the MVP model,
$(z_{t}^{(j)},z_{t}^{(k)})\mid\theta_{jk},y,X$ has a bivariate Gaussian
distribution with mean
$\mu_{t}=(x_{t}^{\mathrm{\scriptscriptstyle{T}}}\beta_{j},x_{t}^{\mathrm{\scriptscriptstyle{T}}}\beta_{k})$,
unit variance and correlation $\sigma_{jk}$. Since the posterior distributions
of the regression coefficients are Gaussian, we can easily marginalize over
them to obtain that
$(z_{t}^{(j)},z_{t}^{(k)})\mid\sigma_{jk}\sim{\mathrm{N}}(\tilde{\mu}_{t},\tilde{\Sigma}_{t})$,
where
$\tilde{\mu}_{t}=(x_{t}^{\mathrm{\scriptscriptstyle{T}}}\hat{\beta}_{j},x_{t}^{\mathrm{\scriptscriptstyle{T}}}\hat{\beta}_{k})$
and the diagonal elements of $\tilde{\Sigma}_{t}$ are
$1+x_{t}^{\mathrm{\scriptscriptstyle{T}}}H_{j}x_{t}$,
$1+x_{t}^{\mathrm{\scriptscriptstyle{T}}}H_{k}x_{t}$. To obtain the
distribution of $(z_{t}^{(j)},z_{t}^{(k)})\mid y,X$, it remains to marginalize
over $\sigma_{jk}$ where
$\sigma_{jk}\sim{\mathrm{N}}(\hat{\sigma}_{jk},s_{jk}^{2})$. In Section 3.3,
we show that $\Pi(\sigma_{jk}\mid
y,X)\overset{d}{\to}\delta_{\sigma_{jk}^{*}}$; a consequence of the total
variation convergence of Theorem 3.7 where $\sigma_{jk}^{*}$ is the true value
of the correlation between outcomes $j$ and $k$. Hence,
$\hat{\sigma}_{jk}\overset{P}{\to}\sigma_{jk}^{*}$. Write the marginal density
$(z_{t}^{(j)},z_{t}^{(k)})\mid\sigma_{jk},y,X$ as the expectation over
$\Pi(\sigma_{jk}\mid y,X)$, that is $p(z_{t}^{(j)},z_{t}^{(k)})\mid
y,X)=E_{\sigma_{jk}}\,p(z_{t}^{(j)},z_{t}^{(k)})\mid\sigma_{jk},y,X)$. Then by
Assumption 3.2 the conditional density
$p(z_{t}^{(j)},z_{t}^{(k)})\mid\sigma_{jk},y,X)$ is bounded and integrable.
Thus, by the dominated convergence theorem,
$p(z_{t}^{(j)},z_{t}^{(k)})\mid\sigma_{jk},y,X)\to
p(z_{t}^{(j)},z_{t}^{(k)})\mid\sigma_{jk}^{*},y,X)$ of which a reasonable
approximation is provided by
$p(z_{t}^{(j)},z_{t}^{(k)})\mid\hat{\sigma}_{jk},y,X)$. With all the above
ingredients, we approximate $p(z_{t}^{(j)},z_{t}^{(k)}\mid y,X)$ by
${\mathrm{N}}(\tilde{\mu}_{t},\hat{\Sigma}_{t})$, where $\tilde{\mu}_{t}$ is
as defined before and $\hat{\Sigma}_{t}$ has the same diagonal elements as
$\tilde{\Sigma}_{t}$ and the off-diagonal elements are $\hat{\sigma}_{jk}$.
Finally, a predictive sample for $(y_{t}^{(j)},y_{t}^{(k)})$ can be drawn by
first drawing the latent variables
$(z_{t}^{(j)},z_{t}^{(k)})\sim{\mathrm{N}}(\tilde{\mu}_{t},\hat{\Sigma}_{t})$
and then thresholding these latent variables. In our numerical experiments,
this provided reliable approximations to the predictive distribution. For
example, when $(n,p,q)=(150,5,2)$ and with 50 test points, the approximation
maintained an average 1% error rate in estimating the posterior predictive
mean for different values of the underlying correlation coefficient in
$[-1,1]$; here we treated the posterior predictive mean computed using the
data augmented MCMC sampler as our benchmark.

##  S.7 Runtime

We compare the runtime of the proposed method with the recently developed fMVP
method (Pichler and Hartig, 2020) which uses the algorithm from Chen et al.
(2018) to compute Gaussian orthant probabilities in parallel and imposes an
elastic net penalty on the correlation matrix for handling large numbers of
outcomes. We leave out TSF of Ting et al. (2022) as their code is not publicly
available (authors declined our request) and the code used to produce results
for their method in the previous subsection is not optimized. In our runtime
analysis we vary the sample size $n$ over a grid from 50 to 500 with
increments of 50 and for each sample size we consider 3 different choices of
the number of outcomes: $q=0.1n,0.3n,0.5n$. The number of covariates $p=5$ and
the true regression matrix $B^{*}$ is simulated according the “dense” case
listed in the previous subsection. The true correlation matrix is set as
$\Sigma^{*}=\rho 11^{\mathrm{\scriptscriptstyle{T}}}+(1-\rho)\mathrm{I}_{q}$
with $\rho=0.5$. For each value of $(n,q)$ we repeated the procedures 20 times
and report the average runtime and average error where the error for each run
is computed as
$\|\hat{B}-B^{*}\|_{F}/pq+\|\hat{\Sigma}-\Sigma^{*}\|_{F}/q^{2}$. The runtime
results are reported in Figure S5, based on analyses conducted on an Intel
i7-8700K CPU with 3.7 GHz processor; the corresponding plot showing the errors
are provided in the supplement. Evidently, bigMVP improves upon the computing
time by an order of magnitude while also performing better in estimation
accuracy, especially in small sample sizes. In addition, unlike the competing
fMVP method, bigMVP provides standard errors and uncertainty intervals
(credible intervals in our case).

Figure S5: Runtime of the hierarchical (red solid), non-hierarchical (green
dashed) versus the fast regularized method (blue dotted) from Pichler and
Hartig (2020) for different numbers of data points $n$ and binary outcomes
$q$. Figure S6: Estimation error of the hierarchical (red solid), non-
hierarchical (green dashed) versus the fast regularized method (blue dotted)
from Pichler and Hartig (2020) for different choices of the number of data
points $n$ and the number of binary outcomes $q$.

##  S.8 Simulation tables and plots

In this section we report results from additional simulation experiments when
the model is misspecified. Specifically, in Figures S7 and S8, we plot
$\|B^{*}-\hat{B}\|_{F}/pq$ (E1) and $\|\Sigma^{*}-\hat{\Sigma}\|_{F}/q^{2}$
(E2) when data are generated from
$t_{10}({B^{*}}^{\mathrm{\scriptscriptstyle{T}}}x,\Sigma^{*})$ and a MVP model
is fitted to the data. In Figure S9, the length of 95% credible intervals are
shown when $n=200$ and $q=10,15,20$, and data are generated under an MVP
model. Tables S.1 and S.2 contains numerical results corresponding to Figures
2 and 3 in the main document.

Figure S7: Logarithm of Frobenius errors in estimating the matrix of
regression coefficients $B^{*}$ when the sample size $n=200$, number of
covariates $p=5$ and the number of binary responses considered are
$q=10,15,20$. The latent variables $z\sim
t_{10}({B^{*}}^{\mathrm{\scriptscriptstyle{T}}}x,\Sigma^{*})$. Figure S8:
Logarithm of Frobenius errors in estimating the matrix of correlation
coefficients $B^{*}$ when the sample size $n=200$, number of covariates $p=5$
and the number of binary responses considered are $q=10,15,20$. The latent
variables $z\sim
t_{10}({B^{*}}^{\mathrm{\scriptscriptstyle{T}}}x,\Sigma^{*})$.  Figure S9:
Comparison of widths of 95% credible/confidence intervals for the correlation
coefficients of bigMVP, bigMVPh and TSF for $n=200$

|  | (Dense, Factor) | (Rare, Factor) | (Dense, Block) | (Rare, Block)
---|---|---|---|---|---
|  | bigMVP | bigMVPh | TSF | bigMVP | bigMVPh | TSF | bigMVP | bigMVPh | TSF | bigMVP | bigMVPh | TSF
$n=200$ | $q=10$ | 92.11 | 93.31 | 88.04 | 98.77 | 98.87 | 35.24 | 93.41 | 93.91 | 81.56 | 96.53 | 96.02 | 58.04
$q=15$ | 93.3 | 93.18 | 87.32 | 97.44 | 97.79 | 47.82 | 93.16 | 93.2 | 87.19 | 97 | 97.42 | 49.51
$q=20$ | 93.89 | 93.81 | 87.63 | 98 | 97.67 | 43.34 | 92.95 | 92.95 | 87.84 | 98.27 | 98.05 | 42.67
$n=500$ | $q=10$ | 94.28 | 94.27 | 92.93 | 95.06 | 95.08 | 71.27 | 94.37 | 94.49 | 91.87 | 95.37 | 95.28 | 71.53
$q=15$ | 94.09 | 94.2 | 92.39 | 95.98 | 95.89 | 65.67 | 95.18 | 95.05 | 92.93 | 95.18 | 95.05 | 70.03
$q=20$ | 93.75 | 93.78 | 91.17 | 95.67 | 95.57 | 66.45 | 93.76 | 93.87 | 92.37 | 94.86 | 94.72 | 69.18

Table S.1: Comparison of bigMVP and bigMVPh versus TSF in terms coverage
probabilities. Here for each case, average coverage probabilities are reported
across all $q(q-1)/2$ correlation coefficients.

|  |  | (Dense, Factor) | (Rare, Factor) | (Dense, Block) | (Rare, Block)
---|---|---|---|---|---|---
|  |  | bigMVP | bigMVPh | HMC | TSF | VB | fMVP | bigMVP | bigMVPh | HMC | TSF | VB | fMVP | bigMVP | bigMVPh | HMC | TSF | VB | fMVP | bigMVP | bigMVPh | HMC | TSF | VB | fMVP
$n=200$ | q=10 | E1 | 0.025 | 0.025 | 0.035 | 0.033 | 0.074 | 0.037 | 0.034 | 0.034 | 0.77 | 0.080 | 0.209 | 0.215 | 0.023 | 0.022 | 0.032 | 0.030 | 0.123 | 0.129 | 0.034 | 0.034 | 1.10 | 0.071 | 0.246 | 0.21
E2 | 0.023 | 0.023 | 0.022 | 0.027 | 0.028 | 0.024 | 0.034 | 0.034 | 0.032 | 0.072 | 0.034 | 0.155 | 0.023 | 0.022 | 0.018 | 0.027 | 0.03 | 0.155 | 0.030 | 0.031 | 0.025 | 0.071 | 0.048 | 0.159
$q=15$ | E1 | 0.019 | 0.019 | 0.026 | 0.025 | 0.067 | 0.023 | 0.026 | 0.026 | 0.724 | 0.064 | 0.107 | 0.057 | 0.019 | 0.019 | 0.037 | 0.025 | 0.07 | 0.027 | 0.026 | 0.026 | 1.77 | 0.062 | 0.131 | 0.055
E2 | 0.015 | 0.015 | 0.014 | 0.018 | 0.023 | 0.015 | 0.022 | 0.022 | 0.021 | 0.048 | 0.024 | 0.023 | 0.016 | 0.015 | 0.011 | 0.021 | 0.015 | 0.015 | 0.019 | 0.020 | 0.015 | 0.050 | 0.016 | 0.016
$q=20$ | E1 | 0.016 | 0.016 | 0.031 | 0.027 | 0.053 | 0.019 | 0.024 | 0.023 | 0.25 | 0.059 | 0.098 | 0.045 | 0.017 | 0.017 | 0.025 | 0.023 | 0.047 | 0.018 | 0.023 | 0.023 | 1.16 | 0.058 | 0.089 | 0.043
E2 | 0.012 | 0.012 | 0.011 | 0.015 | 0.015 | 0.011 | 0.017 | 0.016 | 0.016 | 0.036 | 0.017 | 0.017 | 0.012 | 0.012 | 0.009 | 0.014 | 0.01 | 0.009 | 0.014 | 0.015 | 0.011 | 0.037 | 0.008 | 0.008
$q=100$ | E1 | 0.007 | 0.007 |  |  | 0.016 | 0.015 | 0.009 | 0.008 |  |  | 0.021 | 0.017 | 0.007 | 0.007 |  |  | 0.013 | 0.011 | 0.009 | 0.009 |  |  | 0.026 | 0.024
E2 | 0.002 | 0.002 |  |  | 0.008 | 0.007 | 0.003 | 0.003 |  |  | 0.008 | 0.008 | 0.002 | 0.002 |  |  | 0.009 | 0.008 | 0.002 | 0.002 |  |  | 0.005 | 0.008
$q=200$ | E1 | 0.005 | 0.005 |  |  | 0.006 | 0.005 | 0.005 | 0.005 |  |  | 0.005 | 0.005 | 0.006 | 0.006 |  |  | 0.006 | 0.005 | 0.005 | 0.005 |  |  | 0.007 | 0.006
| E2 | 0.001 | 0.001 |  |  | 0.005 | 0.004 | 0.001 | 0.001 |  |  | 0.005 | 0.005 | 0.001 | 0.001 |  |  | 0.004 | 0.004 | 0.001 | 0.001 |  |  | 0.006 | 0.004
$n=500$ | q=10 | E1 | 0.016 | 0.016 | 0.018 | 0.017 | 0.035 | 0.020 | 0.021 | 0.021 | 0.047 | 0.039 | 0.082 | 0.069 | 0.015 | 0.014 | 0.016 | 0.016 | 0.072 | 0.027 | 0.021 | 0.021 | 0.033 | 0.030 | 0.137 | 0.071
E2 | 0.014 | 0.013 | 0.013 | 0.016 | 0.032 | 0.025 | 0.028 | 0.028 | 0.026 | 0.050 | 0.050 | 0.050 | 0.015 | 0.014 | 0.013 | 0.016 | 0.028 | 0.024 | 0.028 | 0.029 | 0.021 | 0.051 | 0.029 | 0.029
$q=15$ | E1 | 0.014 | 0.013 | 0.016 | 0.015 | 0.055 | 0.018 | 0.017 | 0.017 | 0.099 | 0.028 | 0.107 | 0.052 | 0.012 | 0.012 | 0.013 | 0.012 | 0.042 | 0.011 | 0.016 | 0.016 | 0.060 | 0.035 | 0.108 | 0.055
E2 | 0.009 | 0.009 | 0.008 | 0.010 | 0.022 | 0.013 | 0.019 | 0.019 | 0.018 | 0.034 | 0.023 | 0.023 | 0.010 | 0.010 | 0.008 | 0.010 | 0.012 | 0.012 | 0.017 | 0.017 | 0.011 | 0.036 | 0.014 | 0.014
$q=20$ | E1 | 0.010 | 0.010 | 0.015 | 0.014 | 0.056 | 0.017 | 0.014 | 0.014 | 0.025 | 0.024 | 0.094 | 0.047 | 0.011 | 0.010 | 0.011 | 0.011 | 0.047 | 0.012 | 0.014 | 0.014 | 0.018 | 0.024 | 0.099 | 0.048
E2 | 0.007 | 0.007 | 0.008 | 0.009 | 0.018 | 0.012 | 0.015 | 0.014 | 0.014 | 0.024 | 0.017 | 0.016 | 0.007 | 0.007 | 0.006 | 0.010 | 0.009 | 0.007 | 0.012 | 0.011 | 0.009 | 0.027 | 0.010 | 0.010
$q=100$ | E1 | 0.005 | 0.004 |  |  | 0.012 | 0.009 | 0.005 | 0.005 |  |  | 0.016 | 0.012 | 0.004 | 0.004 |  |  | 0.010 | 0.009 | 0.005 | 0.005 |  |  | 0.019 | 0.018
E2 | 0.001 | 0.001 |  |  | 0.006 | 0.005 | 0.003 | 0.003 |  |  | 0.007 | 0.007 | 0.002 | 0.001 |  |  | 0.004 | 0.005 | 0.003 | 0.002 |  |  | 0.006 | 0.006
$q=200$ | E1 | 0.003 | 0.003 |  |  | 0.005 | 0.005 | 0.003 | 0.003 |  |  | 0.011 | 0.007 | 0.004 | 0.003 |  |  | 0.009 | 0.008 | 0.003 | 0.003 |  |  | 0.010 | 0.008
| E2 | 0.0007 | 0.0007 |  |  | 0.004 | 0.004 | 0.001 | 0.001 |  |  | 0.005 | 0.003 | 0.0008 | 0.0008 |  |  | 0.006 | 0.006 | 0.001 | 0.001 |  |  | 0.005 | 0.005

Table S.2: Comparison of bigMVP, bigMVPh versus TSF, HMC, VB and fMVP in terms
of error in estimating the coefficient matrix $B^{*}$ (E1) and correlation
matrix $\Sigma^{*}$ (E1). In particular,
$\mathrm{E}1=\|B^{*}-\hat{B}\|_{F}/pq$ and
$\mathrm{E}2=\|\Sigma^{*}-\hat{\Sigma}\|_{F}/q^{2}$.

##  S.9 Auxiliary results

######  Lemma S.9.1 (Lemma S6.2 of Miller (2021)).

Let $E\subseteq\mathbb{R}^{D}$ be open and convex, and let
$f_{n}:E\rightarrow\mathbb{R}$ for $n\in\mathbb{N}$. For any $k\in\mathbb{N}$,
if each $f_{n}$ has continuous $k$th-order derivatives and
$\left(f_{n}^{(k)}\right)$ is uniformly bounded, then
$\left(f_{n}^{(k-1)}\right)$ is equi-Lipschitz.

######  Lemma S.9.2.

Let $\\{x_{i}\\}$ be a sequence of $p$-dimensional vectors where each
$x_{i}\in\mathbb{R}^{p}$ and $\lim_{n\rightarrow\infty}n^{-1}\Sigma_{i=1}^{n}$
$x_{i}x_{i}^{\prime}$ is a finite nonsingular matrix. For any $p\times p$
square matrix $A$, we let its trace $\operatorname{tr}(A)$ be the sum of its
diagonal elements, i.e., $tr(A)=\sum_{i=1}^{p}A_{ii}$. Let
$\\{a_{i}\\},\\{b_{i}\\}$ be real number sequences and $\\{b_{i}\\}$ are all
positive. Then

$\sum_{i=1}^{n}(a_{i}x_{i})^{T}\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1}\sum_{i=1}^{n}(a_{i}x_{i})\leq
p\max_{i}\frac{a_{i}^{2}}{b_{i}}$

###### Proof.

We have,

$\displaystyle\frac{1}{n}\sum_{i=1}^{n}(a_{i}x_{i})^{T}\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1}\sum_{i=1}^{n}(a_{i}x_{i})$
$\displaystyle=$
$\displaystyle\frac{1}{n}\operatorname{tr}\left[\sum_{i=1}^{n}(a_{i}x_{i})^{T}\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1}\sum_{i=1}^{n}(a_{i}x_{i})\right]$
$\displaystyle=$
$\displaystyle\frac{1}{n}\operatorname{tr}\left[\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1}\sum_{i=1}^{n}(a_{i}x_{i})\sum_{i=1}^{n}(a_{i}x_{i})^{T}\right]$
$\displaystyle=\frac{1}{n}\operatorname{tr}\left[\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1/2}\sum_{i=1}^{n}(a_{i}x_{i})\sum_{i=1}^{n}(a_{i}x_{i})^{T}\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1/2}\right]$

Also, the matrix

$\displaystyle
A=\sum_{i=1}^{n}\left[a_{i}x_{i}-\sum_{i=1}^{n}\frac{(a_{i}x_{i})}{n}\right]\left[a_{i}x_{i}-\sum_{i=1}^{n}\frac{(a_{i}x_{i})}{n}\right]^{T}$
$\displaystyle=\sum_{i=1}^{n}a_{i}^{2}x_{i}x_{i}^{T}-\frac{1}{n}\sum_{i=1}^{n}(a_{i}x_{i})\sum_{i=1}^{n}(a_{i}x_{i})^{T}$

is positive definite. Plugging the preceding display back to (LABEL:eq:trace),
we have

$\displaystyle\frac{1}{n}\operatorname{tr}\left[\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1/2}\sum_{i=1}^{n}(a_{i}x_{i})\sum_{i=1}^{n}(a_{i}x_{i})^{T}\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1/2}\right]$
$\displaystyle\leq\operatorname{tr}\left[\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1/2}\sum_{i=1}^{n}a_{i}^{2}x_{i}x_{i}^{T}\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1/2}\right]$
$\displaystyle=\operatorname{tr}\left[\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1}\sum_{i=1}^{n}a_{i}^{2}x_{i}x_{i}^{T}\right]$
$\displaystyle\leq\max_{i}\frac{a_{i}^{2}}{b_{i}}\operatorname{tr}\left[\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)^{-1}\left(\sum_{i=1}^{n}b_{i}x_{i}x_{i}^{T}\right)\right]$
$\displaystyle=p\max_{i}\frac{a_{i}^{2}}{b_{i}}$

∎

##  S.10 Some Useful Quantities

Here, we present the detailed expressions of quantities used in Section 3. Set
$\widehat{\Sigma}_{12}=\\{(1,r_{i1}r_{i2}\sigma^{*}_{12})^{\mathrm{\scriptscriptstyle{T}}};(r_{i1}r_{i2}\sigma^{*}_{12},1)^{\mathrm{\scriptscriptstyle{T}}}\\}$
and $\mu_{1}^{*}=x_{i}^{\mathrm{\scriptscriptstyle{T}}}\beta_{1}^{*}$,
$\mu_{2}^{*}=x_{i}^{\mathrm{\scriptscriptstyle{T}}}\beta_{2}^{*}$. We have,

$R^{12}_{1}=\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}\mathrm{Cov}_{\theta^{*}}\left(\frac{\partial\ell_{i}^{1}}{\partial\beta_{1}},\frac{\partial\ell^{12}_{i}}{\partial\sigma_{12}}\right)=0$
(S.14)

where

$\displaystyle\mathrm{Cov}_{\theta^{*}}\left(\frac{\partial\ell_{i}^{1}}{\partial\beta_{1}},\frac{\partial\ell^{12}_{i}}{\partial\sigma_{12}}\right)$
$\displaystyle=E_{\theta^{*}}\left(\frac{\partial\ell_{i}^{1}}{\partial\beta_{1}}\frac{\partial\ell^{12}_{i}}{\partial\sigma_{12}}\right)$
$\displaystyle=E_{\theta^{*}}\left[\frac{\phi(r_{i1}\mu^{*}_{1})}{\Phi(r_{i1}\mu^{*}_{1})}r_{i1}\times\frac{\phi_{\widehat{\Sigma}_{12}}}{\Phi_{\widehat{\Sigma}_{12}}}r_{i1}r_{i2}x_{i}\right]$
$\displaystyle=E_{\theta^{*}}\left[\frac{\phi(r_{i1}\mu^{*}_{1})}{\Phi(r_{i1}\mu^{*}_{1})}\frac{\phi_{\widehat{\Sigma}_{12}}}{\Phi_{\widehat{\Sigma}_{12}}}r_{i2}x_{i}\right]$
$\displaystyle=\frac{\phi(\mu^{*}_{1})}{\Phi(\mu^{*}_{1})}\phi_{[1,\sigma_{12}^{*};\sigma_{12}^{*},1]}(\mu^{*}_{1},\mu^{*}_{2})x_{i}-\frac{\phi(\mu_{1}^{*})}{\Phi(\mu_{1}^{*})}\phi_{[1,-\sigma_{12}^{*};-\sigma_{12}^{*},1]}(\mu_{1}^{*},-\mu_{2}^{*})x_{i}$
$\displaystyle+\frac{\phi(-\mu_{1}^{*})}{\Phi(-\mu_{1}^{*})}\phi_{[1,-\sigma_{12}^{*};-\sigma_{12}^{*},1]}(-\mu_{1}^{*},\mu_{2}^{*})x_{i}-\frac{\phi(-\mu_{1}^{*})}{\Phi(-\mu_{1}^{*})}\phi_{[1,\sigma_{12}^{*};\sigma_{12}^{*},1]}(-\mu_{1}^{*},-\mu_{2}^{*})x_{i}$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}0$

The equation (a) in the above display is because of
$\phi_{[1,\sigma_{12};\sigma_{12},1]}(\mu_{1},\mu_{2})=\phi_{[1,-\sigma_{12};-\sigma_{12},1]}(\mu_{1},-\mu_{2})$.
Similarly,
$R^{12}_{2}=\lim_{n\rightarrow\infty}\frac{1}{n}\mathrm{Cov}_{\theta^{*}}\left(\frac{\partial\ell^{2}_{i}}{\partial\beta_{2}},\frac{\partial\ell^{12}_{i}}{\partial\sigma_{12}}\right)=0$.
For $V^{12}$, we have

$V^{12}=\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}\mathrm{Cov}\left[\frac{\phi(r_{i1}\mu_{1}^{*})}{\Phi(r_{i1}\mu_{1}^{*})}r_{i1},\frac{\phi(r_{i2}\mu_{2}^{*})}{\Phi(r_{i2}\mu_{2}^{*})}r_{i2}\right]$
(S.15)

where

$\displaystyle\mathrm{Cov}\left[\frac{\phi(r_{i1}\mu_{1}^{*})}{\Phi(r_{i1}\mu_{1}^{*})}r_{i1},\frac{\phi(r_{i2}\mu_{2}^{*})}{\Phi(r_{i2}\mu_{2}^{*})}r_{i2}\right]=E\left[\frac{\phi(r_{i1}\mu_{1}^{*})}{\Phi(r_{i1}\mu_{1}^{*})}\frac{\phi(r_{i2}\mu_{2}^{*})}{\Phi(r_{i2}\mu_{2}^{*})}r_{i1}r_{i2}\right]$
$\displaystyle=\sum_{r_{i1}=-1,1}\sum_{r_{i2}=-1,1}r_{i1}r_{i2}\frac{\phi(r_{i1}x_{i}^{T}\beta^{*}_{1})}{\Phi(r_{i1}x_{i}^{T}\beta^{*}_{1})}\frac{\phi(r_{i2}x_{i}^{T}\beta^{*}_{2})}{\Phi(r_{i2}x_{i}^{T}\beta^{*}_{2})}\Phi_{\widehat{\Sigma}_{12}}(r_{i1}x_{i}^{T}\beta^{*}_{1},r_{i2}x_{i}^{T}\beta^{*}_{2}),$

with
$\widehat{\Sigma}_{12}=\\{(1,r_{i1}r_{i2}\sigma^{*}_{12})^{\mathrm{\scriptscriptstyle{T}}};(r_{i1}r_{i2}\sigma^{*}_{12},1)^{\mathrm{\scriptscriptstyle{T}}}\\}$.
Now we turn to $Q^{12}_{1}$.We have

$Q^{12}_{1}=\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}-E_{\theta^{*}}\left(\frac{\partial^{2}{\ell^{12}_{i}}}{\partial{\sigma_{12}}\partial{\beta_{1}^{T}}}\right)$

where

$\displaystyle-
E_{\theta^{*}}\left(\frac{\partial^{2}{\ell^{12}_{i}}}{\partial{\sigma_{12}}\partial{\beta_{1}^{T}}}\right)=\mathrm{Cov}_{\theta^{*}}\left(\frac{\partial\ell^{12}_{i}}{\partial\sigma_{12}},\frac{\partial\ell^{12}_{i}}{\partial\beta_{1}}\right)$
$\displaystyle=\mathrm{Cov}\left(\frac{\phi_{\widehat{\Sigma}_{12}}}{\Phi_{\widehat{\Sigma}_{12}}}r_{i1}r_{i2},\frac{\phi(r_{i1}\mu_{1}^{*})\Phi\left(r_{i2}\frac{\mu_{2}^{*}-\sigma_{12}^{*}\mu_{1}^{*}}{\sqrt{1-(\sigma_{12}^{*})^{2}}}\right)}{\Phi_{\widehat{\Sigma}_{12}}}r_{i1}x_{i}^{T}\right)$
$\displaystyle=\sum_{r_{i1}=-1,1}\sum_{r_{i2}=-1,1}r_{i2}\phi(r_{i1}x_{i}^{T}\beta^{*}_{1})\Phi\left(r_{i2}\frac{x_{i}^{T}\beta^{*}_{2}-\sigma_{12}^{*}x_{i}^{T}\beta^{*}_{1}}{\sqrt{1-(\sigma_{12}^{*})^{2}}}\right)\frac{\phi_{\widehat{\Sigma}_{12}}(r_{i1}x_{i}^{T}\beta^{*}_{1},r_{i2}x_{i}^{T}\beta^{*}_{2})}{\Phi_{\widehat{\Sigma}_{12}}(r_{i1}x_{i}^{T}\beta^{*}_{1},r_{i2}x_{i}^{T}\beta^{*}_{2})}x_{i}^{T}$
(S.16)

Similarly,

$\displaystyle Q^{12}_{2}=$
$\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{i=1}^{n}\sum_{r_{i1}=-1,1}\sum_{r_{i2}=-1,1}r_{i1}\phi(r_{i2}x_{i}^{T}\beta^{*}_{2})\Phi\left(r_{i1}\frac{x_{i}^{T}\beta^{*}_{1}-\sigma_{12}^{*}x_{i}^{T}\beta^{*}_{2}}{\sqrt{1-(\sigma_{12}^{*})^{2}}}\right)\frac{\phi_{\widehat{\Sigma}_{12}}(r_{i1}x_{i}^{T}\beta^{*}_{1},r_{i2}x_{i}^{T}\beta^{*}_{2})}{\Phi_{\widehat{\Sigma}_{12}}(r_{i1}x_{i}^{T}\beta^{*}_{1},r_{i2}x_{i}^{T}\beta^{*}_{2})}x_{i}^{T}$
(S.17)


††thanks: $\star$ _equal first authorship_ 1 Massachusetts Institute of
Technology
2 Weizmann Institute of Science 3 Harvard University
4 Fujitsu Research of America, Inc. 5 Yale University 6 Fujitsu Limited
7 Center for Brains, Minds and Machines

# Emergent Neural Network Mechanisms for Generalization to Objects in Novel
Orientations

Avi Cooper$\star$,1,4,5,7, Xavier Boix$\star$,1,4,7, Daniel Harari1,2,7,
Spandan Madan3,7,
Hanspeter Pfister3, Tomotake Sasaki6,7, Pawan Sinha1

###### Abstract

The capability of Deep Neural Networks (DNNs) to recognize objects in
orientations outside the distribution of the training data is not well
understood. We present evidence that DNNs are capable of generalizing to
objects in novel orientations by disseminating orientation-invariance obtained
from familiar objects seen from many viewpoints. This capability strengthens
when training the DNN with an increasing number of familiar objects, but only
in orientations that involve 2D rotations of familiar orientations. We show
that this dissemination is achieved via neurons tuned to common features
between familiar and unfamiliar objects. These results implicate brain-like
neural mechanisms for generalization.

a

b

Figure 1: Learning paradigm and network’s per-orientation accuracy. (a) An
intuitive presentation of our experimental paradigm. The network is trained
with images of certain airplanes at all orientations, constituting a ‘full
experience.’ The network is also trained with images of a small subset of
orientations for other airplanes, constituting a ‘partial experience.’ The
out-of-distribution, or OoD, generalization capacities of the network are
evaluated by measuring the classification accuracy for partially-seen
airplanes at unseen orientations. Our results suggest that OoD generalization
is facilitated by the dissemination of orientation invariance developed for
all orientations for the fully-seen airplanes to the OoD orientations of the
partially-seen airplanes. (b) Left: The learning paradigm employed in this
work. Each column is a sample object instance (here from the airplane dataset)
and each row is a sample orientation. The training set includes all
orientations for fully-seen instances, and a partial set of orientations
(outlined in red) for partially-seen instances (in this example, with the
airplanes’ nose pointing down). The orientations included in the training set
are referred to as in-distribution orientations (pink shading). Orientations
of the partially-seen instances that are not included in the training set are
referred to as OoD (yellow shading). Right: A visualization of per-
orientation-analysis. The set of all orientations can be portrayed as a
square, where proximity relationships between orientations are captured in
their arrangement in the square. (Further details are provided in Fig.2a.) The
in-distribution and OoD sets for partially-seen airplanes are represented in
their respective shadings.

##  1 Introduction

Recognizing objects in novel orientations lies at the heart of biological and
artificial intelligence, as it is a fundamental capacity to understand the
visual world [1, 2]. However, it remains as an outstanding question.

In the realm of artificial systems, Deep Neural Networks (DNNs) have recently
made large strides in learning to recognize objects. However, DNNs’
generalization is often limited to the distribution of images used for
training, known as in-distribution data, and DNNs’ performance tends to
deteriorate when confronted with out-of-distribution (OoD) data. Previous
studies have shown that DNNs perform poorly when objects are presented in
novel orientations, even when learning from large datasets with millions of
examples [3, 4, 5].

An approach to understand the capabilities of DNNs is to leverage the
knowledge gained from studying biological intelligence [6, 7]. In natural
settings, biological intelligent agents observe instances of object categories
from diverse orientations. When encountering a new object instance, these
agents often demonstrate the capacity to accurately identify the object in
different orientations by drawing upon past experiences with similar
instances. There has been a extensive investigation dedicated to the
exploration of human and mammalian perceptual capabilities in the domain of
object recognition in unfamiliar orientations and the neural mechanisms
underlying these cognitive abilities. Studies have revealed that recognition
accuracy varies across novel orientations, with some orientations exhibiting
superior generalization compared to others [8]. Additionally, compelling
evidence suggests that neurons responds to their own specific set of object
features when these are present in the visual field [9, 10]. This neural
tuning has been reported to be invariant to a certain degree from the object’s
orientation [11]. Theoretical frameworks have proposed that such neural
invariance to object orientation forms the basis for the ability to recognize
objects in novel orientations within biological systems [12].

We study DNNs under conditions akin to the operating regime of biological
brains, in which some instances of an object category (_e.g.,_ a ‘Boeing 777
airliner’ is an instance of the ‘airplane’ category) are seen from all
orientations during training (fully-seen instances), while other instances are
only seen in a subset of all orientations (partially-seen instances). During
test time, we evaluate the generalization performance of the networks by
measuring instance classification performance on OoD orientations (_i.e.,_
those orientations not included in the training set) of partially-seen
instances. This is a simple paradigm, inspired by the paradigm by [13], that
facilitates analyzing the impact of several key factors that may influence OoD
generalization, such as the number of fully-seen instances and the in-
distribution orientations of the partially-seen instances. This paradigm
allows us to more precisely characterize performance challenges of DNNs for
OoD orientations. Figure 1 summarizes the paradigm that we follow in this
work.

In this paper, we employ the same analytical tools utilized in the study of
biological brains, with the aim of understanding the DNNs’ generalization
abilities in OoD orientations. Specifically, we investigate questions that
have already been investigated within the realm of biological systems, such as
whether DNNs exhibit different recognition failures across OoD orientations,
and explore whether individual neurons within DNNs display feature tuning for
object recognition while also exhibiting invariance to different orientations.
We aim at explaining the generalization abilities observed in DNNs.

##  2 Results

a | b |  |
---|---|---|---
|  |  |
| 20 Fully Seen | 40 Fully Seen |

c |  |
---|---|---
|  |
$\gamma\approx-\pi$ | $\gamma\approx\frac{-1}{2}\pi$ | $\gamma\approx 0$

Figure 2: Observed generalization patterns in per-orientation accuracy
heatmaps. When trained with a combination of fully-seen instances and
partially-seen instances, DNNs demonstrate the ability to generalize outside
of their training distribution. Generalization behaviors are demonstrated
measuring per-orientation accuracy. (a) All orientations can be described by
three Euler axes ($\alpha,\beta,\gamma$,) and rotations are periodic around
these axes. These properties allow for the visualization of all possible
orientations with an orientation cube, shown here. The orientations contained
within the colored rectangular prism are those orientations of the partially-
seen instances included in training (_i.e.,_ are in-distribution). The in-
distribution orientations differ depending on the experiment. All other
orientations are OoD. (b) Increased network generalization for OoD
orientations, with increased instance diversity (_i.e.,_ number of fully-
seen.) Each cell in the heatmap is the average classification accuracy of the
network for a given value of $\beta$ and $\gamma$, across all values of
$\alpha$. Chance level is 0.02 (2%). (c) Different in-distribution parameters
affect the generalization behaviors. The generalization patterns for a
different span of in-distribution orientations ($-0.25\leq\alpha\leq
0.1,-0.1\leq\beta\leq 0.25,-\pi\leq\gamma<\pi$) as outlined by the purple box.
In this case, each cell is of a given value for $\alpha,\beta,\gamma$.

a
b

Figure 3: Modeling generalization patterns for OoD orientations. The bar plots
show several trends related to DNN OoD classification patterns. The trends are
measured under the various controls, including in-distribution orientations
conditions ($\alpha,\gamma,\beta,\alpha^{\prime}$) and object category, which
is either a single object as in Airplane, SM, Car, or transfer across two
categories, when the fully-seen instances are of a different category than the
partially-seen instances as in Airplane $\rightarrow$ SM and vice versa. These
transfer cases are visually separated from the other cases. (a) Network
generalization for OoD orientations increases with increasing number of fully
seen (blue shading.) This trend holds across object category and in-
distribution orientations conditions. (b) Top: We introduce a predictive model
for OoD orientation generalization (black — “All Components”) which is highly
predictive of experimental results, with greater than 0.8 Pearson Correlation
Coefficient for all experimental controls. (Results are shown for experiments
with 40 fully-seen instances.) Null hypothesis predictive models, including
“Random Uniform” and “In-Distribution,” have very low correlation
coefficients. We also ablate our predictive model, including only some sub-
components, like only-“Small Angle”, only-“In-Plane” or only-“Small Angle +
In-Plane.” These ablated models have lower correlation coefficients than “All
Components,” and vary in relation to one another depending on the experimental
condition. Bottom: We isolate the predictive power of the only-“In-Plane”
component for all experiments with a range of number of fully-seen. The
increasing predictive power of the “In-Plane” component correlates with
increasing OoD accuracy as the number of fully-seen instances increases. This
suggests that generalization to “In-Plane” orientations drives OoD accuracy.

a
b c

Figure 4: Neuronal analysis, Invariance and Dissemination. (a) An intuitive
visualization of neural activity. Each square is the response of a single
neuron to the airplane instance that most highly activates it portrayed in two
ways: 1) the top-8 images that most highly activate the neuron (in no
particular order), 2) the heatmap of the per-orientation normalized neural
activity for the airplane instance. Neurons tend to exhibit patterns of
activation related to the patterns of generalization behavior (Figs. 2b for
example,) and are invariant to a range of orientations that respect the
partitioning of OoD orientations. The patterns of activation are similar
between the fully-seen instance that most highly activates the neuron and the
partially-seen instance that most highly activates it due to shared visual,
part and semantic features between these instances. Several randomly sampled
penultimate layer neurons, arranged into columns, demonstrate that these
findings apply to many neurons. (b) Averaging the activations in the
partitioned regions (‘seed’, generalizable, non-generalizable) and computing
the invariances (defined here: Eq.6) between ‘seed’ and OoD regions captures
overall generalization in the network. Plotting the generalization metrics
against accuracy for those regions demonstrates a clear correlation between
increasing invariance and increasing OoD classification accuracy (_i.e.,_
partially-seen instances in OoD orientations.) (c) Plotting fully-seen
invariance against partially-seen invariance also yields a tight correlation,
suggesting that dissemination of invariance from fully-seen to partially-seen
instances enables increasing generalization in OoD orientations of partially-
seen instances.

#### Overview.

In exploration of the generalization capabilities of DNNs to novel
orientations, our computational experiments show that the OoD orientation
space is divided between generalizable and non-generalizable orientations, in
terms of the networks behavior. We find this division to be governed by a set
of rules which determine a partitioning of the orientation space, given a
‘seed’ of orientations seen by the DNNs at train time. Among several
partitioning rules, we identify a rather intuitive one: small 3D perturbations
around seen orientations will be included in the highly-generalizable
partition. We find other rules to be more surprising, including that the
highly-generalizable partition also consists of shape and silhouette
preserving rotations, such as in-plane (_i.e.,_ 2D) rotations and flips along
axes of symmetry of the seen orientations. Our analysis shows that the
predicted network behavior induced by these partitioning rules is highly
correlated with the measured network behavior under a variety of training
regimes. We further explore the DNNs internal representations and identify
neuronal mechanisms in the network, that allow the dissemination of
orientation-invariance from familiar objects to novel objects and
orientations.

#### Per-orientation accuracy heatmaps.

A detailed inspection of the network’s generalization capability for OoD
orientations is brought forward by introducing per-orientation accuracy
heatmaps. In brief, the continuous space of object orientations, represented
by Euler angles, is discretized into cubelets (Fig. 2a). For each cubelet the
network’s performance is evaluated in terms of the classification accuracy
$\Psi(\bm{\theta})$, where $\bm{\theta}$ is an orientation of interest.
Accuracy heatmaps are 2D projections across a specified dimension of the full
accuracy orientation cube (Fig. 2b,c; Methods). These heatmaps reveal a
reproducible pattern of generalization in the form of increased classification
accuracy for OoD (_i.e.,_ novel) orientations. For example, Figure 2b shows
that for ‘seed’ orientations at the center of the heatmap (red box), the
network (in this experiment - ResNet18 [14]; see Methods) yields the highest
accuracy (brightest cubelets) for adjacent orientations around the ‘seed’,
depicting small 3D perturbations of the ‘seed’ orientations. Further
inspection of the heatmap reveals other orientations, in this example,
brighter cubelets forming the figure ‘8’ (stretching ‘seed’ sideways and along
the heatmap’s boundaries, enclosing two darker ‘holes’), for which the network
performs better than for the rest of the OoD orientations. These orientation
mainly depict in-plane rotations of the ‘seed’ orientations.

When considering the average classification accuracy across all OoD
orientations, our experiments reproduce previous results. In particular, we
can reliably quantify the effect of data diversity on the OoD generalization,
as the amount of training examples is kept constant with our learning
paradigm. Figure 3a clearly shows an increase in OoD accuracy as data
diversity (_i.e.,_ the number of fully-seen instances) increases, under
various conditions, including different ‘seed’ orientations, different image
datasets and across datasets. The accuracy heatmaps provide a complementary
means of assessment to the overall average accuracy measure, depicting the
generalization patterns and indicating which orientations account for the
network increased performance (_e.g.,_ Fig. 2b). The patterns of increased
accuracy depict a partitioning of the orientation space, which reappears for
various ‘seed’ orientations (Fig. 2c), various sizes of the training set and
different object categories (e.g., Airplane, Car, Shepard & Metzler (SM)
objects [15]), as shown in several experiments (Supplementary Figs. S.1, S.2;
Methods).

#### Modeling generalization patterns.

We hypothesize a set of rules which govern the partitioning of the orientation
space into generalizable and non-generalizable orientations. To evaluate this
hypothesis we formulate a model of the partitioning rules, which can be used
to predict the OoD generalization patterns of the network, given a ‘seed’ of
in-distribution orientations. Briefly, the model, denoted by
$f_{\mathbf{w}}(\bm{\theta})$ has three components: $A(\bm{\theta})$, which
captures small angle rotations around $\bm{\theta}$; $E(\bm{\theta})$, which
captures in-plane (2D) rotations; $S(\bm{\theta})$, which captures object
silhouette projections at the orientation $\bm{\theta}$ (see details in
Methods). We evaluate the model’s performance by measuring the Pearson
correlation coefficient $\rho$ between the accuracy of the networks as
measured in our experiments and as predicted by the model,  _i.e.,_
$\rho(\Psi(\bm{\theta}),f_{\mathbf{w}}(\bm{\theta}))$. Figure 3b shows the
predictive power of the model and its components in experiment with different
‘seed’ orientations and several object categories. The model’s component
$A(\bm{\theta})$ (‘small angle’ rotations), is the best predictor for the
network’s OoD behaviour, for highly articulated objects such as the SM
objects. On the other hand, the model’s component $E(\bm{\theta})$ (‘in-plane’
rotations), is a better predictor for non-articulated objects with inherent
symmetries. Further analysis of this component illustrates how generalization
to ‘in-plane’ rotations emerges with the increase in data diversity.

We conducted a large series of experiments under various settings, including
different ‘seed’ orientation distributions, various amounts of training
examples, different object scales, object categories with different levels of
symmetry, image datasets and DNN architectures (see Methods, Supplementary
Figs. S.1, S.2, S.3, S.4, S.5). In all experiments our model highly predicts
the network’s behaviour, indicating that indeed the networks generalization
patterns for OoD orientations follow the model’s partitioning rules. (See S.6
for reports of average accuracy in generalizable and non-generalizable
orientations.) This is true even across categories, when the ‘seed’ is taken
from one category (_e.g.,_ SM) and the ‘fully-seen’ instances are taken from
another (_e.g.,_ Airplane).

#### Individual unit neuronal analysis.

In search of how generalization and dissemination emerge in DNN’s we turn to
analyze the neurons’ activation in the trained networks. We focus on neurons
in the penultimate layer of the network, which are attuned to the highest
level features in the input stimuli, but reflect a consolidated representation
of the entire network for inferring the downstream task (classification in our
simulations).

Figure 4a illustrates activation of individual neurons for stimuli of fully-
seen and partially-seen instances. Each group of images depict input stimuli
of a particular instance for which a particular neuron has the highest
activation, along with the neuron’s per-orientation activation heatmap for the
instance. The patterns seen in the neurons’ activation heatmaps resemble the
partitioning patterns of the accuracy heatmaps shown in Figure 2b. Some
neurons exhibit similar activation patterns for both fully-seen and partially-
seen instances, while others do not.

A quantifiable measure of these neuronal responses can help understanding how
generalization occurs in the network, in particular for the partially-seen
instances observed only at the ‘seed’ orientations during training. Hence, any
generalization to OoD orientations for the partially-seen instances must stem
from the ‘seed’ orientations. We define an activation invariance score in the
range ${[0,1]}$ (Eq. 6) between sets of orientations, in particular between
the ‘seed’ orientations and OoD generalizable orientations or non-
generalizable orientations. The invariance score yields higher values when a
neuron fires for both sets of orientations, and lower values when it fires
only for one set (see details in Methods). We expect that generalization,
reflected by the accuracy level, will correlate with the invariance score.

Figure 4b depicts a scatter plot of the invariance score against the
classification accuracy. Each dot represents an instance at an orientation
set, and the coloring indicates the respective instance set (fully-seen or
partially-seen) and orientation set (generalizable or non-generalizable). As
expected, there is a clear correlation between increasing levels of
classification accuracy and increasing invariance score for the partially-seen
instances. Furthermore, the plot shows a clear partition between generalizable
and non-generalizable orientations with respect to the invariance score, where
significantly higher invariance scores are measured for the generalizable
orientations.

For fully-seen instances (Fig. 4b gray dots), all orientations are in-
distribution, including the ‘seed’, generalizable and non-generalizable.
Therefore, the network easily achieves accuracy at ceiling levels regardless
of the neuronal invariance score. Nevertheless, the fully-seen instances
exhibit the same invariance partitioning between generalizable and non-
generalizable orientations as the partially-seen instances.

Figure 4c depicts a direct comparison between the invariance score of the
fully-seen and partially-seen sets for both the generalizable and non-
generalizable orientations. The partition between generalizable and non-
generalizable orientations is exhibited again — the non-generalizable
invariances are in the bottom left corner, while the generalizable invariances
are in the top right corner. Each point in this plot represents the joint
invariance of fully-seen and partially-seen instances at a given orientation.
The plot clearly shows a tight correlation between the invariance scores of
the fully-seen and partially-seen instances, as most of the points lie within
a band roughly 0.1 units away from the line of parity, $x=y$. This correlation
suggests that an increase in the invariance score of the network at a set of
orientations for the fully-seen instances will be disseminated to the
partially-seen instances (see also S.7).

##  3 Discussion

Our results support the hypothesis that the network disseminates orientation-
invariance of fully-seen instances to partially-seen instances using brain-
like mechanism like those reported by [11, 12]. Neurons are feature detectors
and some of the features that neurons are tuned are features shared between
fully-seen and partially-seen instances (Fig.4a), and therefore the same
neuron helps detecting features for several objects. Features are seen at
multiple orientations in the fully-seen object instances, which enables
neurons to develop orientation invariance. Since features detected by the
neurons are shared among fully-seen and partially-seen instances, invariance
that develop for fully-seen instances are gained “for free” for partially-seen
instances in the same orientations. Note that orientation-invariance is
learned through the fully-seen instances, as increasing the number of fully-
seen instances (while keeping the number of training examples constant)
results in an increase of the DNN capability of disseminating orientation-
invariance obtained from familiar objects. Our results elucidate the intricate
neural processes involved in object recognition and also underscores the
critical role of individual neuron, feature-based representations for OoD
object recognition.

Dissemination of orientation-invariance has been observed for orientations
that appear like 2D rotations (_in-plane_) of in-distribution orientations. In
some cases, when the network relies on the object instance’s silhouette for
recognition, the in-distribution orientations also include orientations that
have the same silhouette as the seen orientations. For non-generalizable
orientations, the network has not developed orientation-invariance with
respect to the seed orientations (demonstrated by the lower invariance score
in our results). It is worth noting that despite the absence of orientation-
invariance, the network is still able to recognize fully-seen instances in
such non-generalizable orientations. This is due to the fact that these
orientations fall within the training distribution and the network has learned
to associate them with their corresponding object instances. However, in the
case of non-generalizable orientations, the dissemination of orientation-
invariance is not feasible. This is even the case when neurons are tuned to
features shared with partially-seen instances, as they do not exhibit
orientation-invariance for these non-generalizable orientations and the
training process does not provide any information to establish associations
with the corresponding object instances. These findings reveal discernible
patterns in the successes and failures of DNNs across diverse orientations.
Such patterns can be effectively characterized and explained through the
analysis of neural activity. This underscores the potential for more
comprehensive analyses of DNNs that transcend the conventional approach of
solely focusing on average accuracy.

A key question arising from our results is to explain why DNNs disseminate
orientation-invariance only to in-plane orientations. All object instances are
distinguishable at all orientations, as evidenced by the high in-distribution
accuracy achieved by the DNNs. Therefore the lack of orientation-invariance
for such non-generalizable orientations is an outcome of the DNN’s learning
process. We speculate that this may be because orientations that are not in-
plane are affected by self-occlusion, which poses a particular challenge for
DNNs. Various efforts have been made to enhance DNNs’ generalization
capabilities to OoD orientations including leveraging preconceived components
for DNNs, such as 3D models of objects [16] or sophisticated sensing
approaches like omnidirectional imaging [17]. However, these approaches rely
on _ad-hoc_ approaches tailored to specific objects and do not address the
fundamental limitations of the DNN learning process in recognizing objects in
OoD orientations. Instead, novel network architectures that extend the
emergent orientation-invariance inherent within networks might allow for
further gains of OoD generalization. Biological agents may overcome the
difficulties associated with recognizing OoD orientations by leveraging the
temporal dimension to associate orientations and learn invariant
representations [18, 19]. The mechanisms that utilize temporal association may
hold fundamental significance, given that they have access to a plentiful
source of training data that does not rely on external guidance and task
specific labels. This data is readily available prior to any visual task and
has the potential to contribute to the emergence of orientation-invariant
representations beyond _in-plane_ orientations.

Previous studies have extensively compared the behavioural and
electrophysiological aspects of brains and DNNs [20]. However, a direct
comparison between these systems alone has limitations in providing insights
into the underlying mechanisms of object recognition in DNNs. This is due to
the possibility that while certain fundamental mechanisms may be shared across
these systems, the manifestation of these fundamental mechanisms can differ at
the behavioral and electrophysiological levels. Our study has provided
compelling evidence of brain-like neural mechanisms in DNNs that facilitate
object recognition in novel orientations, even though these mechanisms are
manifested differently than in biological systems. For instance, while humans
and primates can recognize objects in orientations that are not simply 2D
rotations, this capability is not fully replicated in DNNs. Thus, we can
conclude that the neural mechanisms that have been observed to govern
recognition in biological systems largely apply to DNNs, albeit with distinct
manifestations across these systems. It will be interesting to follow this
line of investigation across biological and artificial systems to envision a
general theory to explain emergent mechanisms in both brains and machines.

##  4 Methods

###  4.1 Per-Orientation Accuracy Visualization

Previous works have typically reported average performance over all
orientations. In contrast, we evaluate the network’s performance for each
orientation across the entire range of orientations. To express an orientation
of an object instance we use $\bm{\theta}:=(\alpha,\beta,\gamma)$, the Euler
angles with respect to the orthogonal axes of a reference coordinate system
${\mathbb{R}}^{3}$ [21], with the convention that $\alpha$ and $\gamma$ are
bounded within $2\pi$ radians, and $\beta$ is bounded within $\pi$ radians. We
define $\Psi(\bm{\theta})\in[0,1]$ to be the network’s average classification
accuracy at an orientation $\bm{\theta}=(\alpha,\beta,\gamma)$ over either the
fully-seen or partially-seen instances.

To facilitate intuition of $\Psi$ we introduce a visual representation of this
function. Since orientations are continuous values and are related spatially
we map the range of bounded values of orientations $(\alpha,\beta,\gamma)$
onto a Cartesian coordinate system, resulting in a cube—the basis of our
visualization. We discretize the continuous space of orientations into
cubelets, which are sub-cubes with a width of $\frac{1}{\\#\textit{Cubelets}}$
of the full range of each respective angle. This approach preserves local
behavior in aggregate analysis. In addition, we outline the range of
orientations which are in-distribution for the partially-seen instances — the
rest are OoD orientations. To illustrate the object orientation at a given
cubelet, we sample one representative image and overlay it onto the heatmap at
the location of the cubelet.

See Fig. 2a, which shows this visual representation scheme, and Figs. 2b and
2c for examples.

###  4.2 Model of DNN Per-Orientation Generalization

In Results we briefly introduce the hypothesis that DNNs are capable of
generalizing to orientations which are small angle rotations of the in-
distribution orientations images and to orientations that are in-plane
relative to the in-distribution images. In this section we formalize this
model.

Recall that we defined $f_{\mathbf{w}}(\bm{\theta})$ as the predictive model
for generalization per each orientation. To measure the goodness of our
prediction, we employ the Pearson correlation coefficient to measure how
closely our model correlates with DNN recognition accuracy,
$\Psi(\bm{\theta})$. We choose this metric because it normalizes data with
respect to amplitude and variance, and therefore measures patterns of behavior
across $\bm{\theta}$ and relative to other $\bm{\theta}$, rather than the
exact performance for every $\bm{\theta}$.

Our model $f_{\mathbf{w}}(\bm{\theta})$ is composed by three components
($A(\bm{\theta})$, $E(\bm{\theta})$ and $S(\bm{\theta})$), which we introduce
next. These three components easily lend themselves to formalization with
Euler’s rotation theorem [21]. The theorem states that any rotation can be
uniquely described by a single axis, represented by a unit vector
$\hat{\mathbf{e}}\in{\mathbb{R}}^{3}$, and an angle of rotation, denoted as
$\phi\in[0,\pi]$ around the axis $\hat{\mathbf{e}}$. We employ this
representation to describe the rotation between an arbitrary orientation of
interest, $\bm{\theta}$, and an orientation in the set of in-distribution,
denoted $\bm{\theta}_{s}\in\Omega_{s}$. We use
$\hat{\mathbf{e}}_{\bm{\theta},\bm{\theta}_{s}}$ and
$\phi_{\bm{\theta},\bm{\theta}_{s}}$ to denote the unit vector (axis) and the
angle of this rotation, respectively.

#### Component 1: Small Angle Rotation, $A(\bm{\theta})$.

The first component of the model captures orientations that are small angle
rotations from the orientations in the training distribution. Visually similar
orientations are those that are arrived at by small rotations from in-
distribution orientations, or small $\phi_{\bm{\theta},\bm{\theta}_{s}}$. We
therefore define the first component $A(\bm{\theta})$ as

$A(\bm{\theta}):=\max_{\bm{\theta}_{s}\in\Omega_{s}}\left|1-\frac{\phi_{\bm{\theta},\bm{\theta}_{s}}}{\pi}\right|\in[0,1].$
(1)

The $\max_{\bm{\theta}_{s}\in\Omega_{s}}$ operator chooses the in-distribution
orientation that is closest to $\bm{\theta}$ of interest.

#### Component 2: In-plane Rotation, $E(\bm{\theta})$.

The second component of the model captures orientations which appear as in-
plane rotations of in-distribution images. Let $\mathbf{c}\in{\mathbb{R}}^{3}$
be the unit vector representing the camera axis. In-plane rotations are those
for which the axis of rotation is parallel to the camera axis. Thus, an
orientation appear as an in-plane rotations of an in-distribution images when
$\mathbf{c}\in{\mathbb{R}}^{3}$ and
$\hat{\mathbf{e}}_{\bm{\theta},\bm{\theta}_{s}}\in{\mathbb{R}}^{3}$ (_i.e.,_
the vector of object instance rotation) are parallel. Taking their standard
inner product yields the proximity to being parallel, which is therefore the
degree to which the rotation is in-plane.

Thus, we define the second component $E(\bm{\theta})$ as follows:

$E(\bm{\theta}):=\max_{\bm{\theta}_{s}\in\Omega_{s}}\left|\mathbf{c}^{\top}\hat{\mathbf{e}}_{\bm{\theta},\bm{\theta}_{s}}\right|\in[0,1],$
(2)

where $\mathbf{c}^{\top}$ denotes the transpose of $\mathbf{c}$.

#### Component 3: Silhouette, $S_{A}(\bm{\theta})$, $S_{E}(\bm{\theta})$.

The third component of the model captures orientations which project object
silhouettes onto the camera that are similar to the silhouettes of the object
when in-distribution — for example, the airplane when viewed from above, and
the silhouette being the airplane viewed from below. These orientations are
defined as a $\pi$ radians rotation around the $\gamma$ axis, which results in
a silhouette orientation. We transform all the in-distribution orientations,
$\Omega_{s}$, in this way, and we call these silhouette in-distribution
orientations $\Omega_{\hat{s}}$. We then compute $S_{A}(\bm{\theta})$ and
$S_{E}(\bm{\theta})$, substituting $\Omega_{\hat{s}}$ for $\Omega_{s}$ in
$A(\bm{\theta})$ and $E(\bm{\theta})$ respectively.

#### Nonlinearities.

The components described above capture a general trend, but do not match the
range of values given by a 0-100% accuracy metric. We therefore fit the
components with a logistic function. The ‘S’-like shape of the logistic
function allows for the highest and lowest values of $E(\bm{\theta})$,
$A(\bm{\theta})$, $S_{A}(\bm{\theta})$ and $S_{E}(\bm{\theta})$ to be close to
the highest and lowest values of $\Psi(\bm{\theta})$. In addition, it allows
for a smooth transition between these highest and lowest values. Most
importantly, the simplicity of the logistic function allows for fitting while
preserving the interpretability of the model components, ensuring that the
models remains related to small angle, in-plane and silhouette rotations. We
employ the following logistic function:

$\sigma(x;(a,b,c))=\frac{1}{1+e^{b(-x^{c}+a)}},$ (3)

where
$x\in\\{E(\bm{\theta}),A(\bm{\theta}),S_{A}(\bm{\theta}),S_{E}(\bm{\theta})\\}$.
$a$ and $b$ translate and scale the values of the predictive components and
$c$ spreads out saturated values of the component.

#### Fitting the Model with Gradient Descent.

The model combines four components $A(\bm{\theta})$, $E(\bm{\theta})$,
$S_{A}(\bm{\theta})$ and $S_{E}(\bm{\theta})$ by taking the sum of their
respective values after applying the logistic function $\sigma$:

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle
f_{\mathbf{w}}(\bm{\theta})$\mbox{}\hfil\\\ :=\hskip
2.84526pt&\sigma(A(\bm{\theta});\mathbf{w}_{A})+\sigma(E(\bm{\theta});\mathbf{w}_{E})+\\\
&\sigma(S_{A}(\bm{\theta});\mathbf{w}_{SA})+\sigma(S_{E}(\bm{\theta});\mathbf{w}_{SE}),\end{split}$
(4)

where $\mathbf{w}$ represents the parameters of the logistic functions
_i.e.,_
$\mathbf{w}=(\mathbf{w}_{A},\mathbf{w}_{E},\mathbf{w}_{SA},\mathbf{w}_{SE})$.
The logistic fitting function is differentiable, and
$f_{\mathbf{w}}(\bm{\theta})$, the linear combination of these logistic
functions, is also differentiable. Further, the Pearson correlation
coefficient is also differentiable. Therefore we employ gradient descent to
fit $\mathbf{w}$ with the Pearson correlation coefficient as the cost
function.

###  4.3 Neural Analysis

In Results we discussed our findings that OoD generalization in the network is
allowed for by dissemination of orientation invariance from fully-seen
instances to partially-seen instances. In this section, we outline the process
by which we quantify several different network invariance metrics. We first
formalize the notation for neural activations for single orientations and for
sets of orientations. We then define the invariance score (Eq. 6). Finally, we
average together many invariance calculations to arrive at the network
invariance metric.

We begin by formalizing our approach to neural activations. In Sec.4.1 we
introduced $\Psi(\bm{\theta})$, the network’s average accuracy at a specific
orientation. We can similarly define the neural activation at a specific
orientation, though we do so with more granularity. Namely, we introduce
$\Phi^{n}_{i}(\bm{\theta})$, which is the average activation of a neuron $n$
from the set of all penultimate-layer neurons $N$ (_i.e.,_ $n\in N$) across
all images of an object instance $i$ from the set of all object instances $I$
(_i.e.,_ $i\in I$) for a given orientation $\bm{\theta}$. We normalize the
activity of each neuron by dividing the activity level of each image by the
maximum activity generated by any image. We exclude any neurons with a maximum
activation of $0$ from further analysis.

Having defined $\Phi$ we note that it is useful to perform analysis not on
single orientations only, but sets of orientations. We demonstrated that under
our experimental conditions, orientations can be partitioned into coherent
subsets — in-distribution and OoD orientations. Further, the OoD orientations
can be partitioned into generalizable orientations, _i.e.,_ those OoD
orientations that the network can generalize to, and non-generalizable
orientations. We refer to the in-distribution, generalizable and non-
generalizable orientation sets as $\mathit{InD}$, $G$ and $\neg G$
respectively. The determination of membership of the generalizable and non-
generalizable orientation sets is as follows: We compute 10% of the maximum
value of $f_{\mathbf{w}}(\bm{\theta})$, the predictive model, in the
experiment with 40 fully-seen instances. All orientations for which $f$ is
greater than the 10% threshold are considered generalizable, otherwise they
are considered non-generalizable. We can now compute the average activation of
a set or orientations. For example, the average activation for a given neuron
$n$ and object instance $i$ of the generalizable orientations is defined in
the following way:

$\bar{\Phi}_{i}^{n}(G)=\frac{1}{|G|}\sum_{\theta\in G}\Phi_{i}^{n}(\theta).$
(5)

The same may be computed for in-distribution and non-generalizable
orientations.

To determine how dissemination occurs in the network, we calculate the degree
of similarity in a neuron’s response to a given instance across different
orientations. Specifically, given a neuron $n$ and instance $i$, we calculate
the similarity between the neuron’s response at an orientation pair
$\Phi^{n}_{i}(\bm{\theta}_{1})$, and $\Phi^{n}_{i}(\bm{\theta}_{2})$, or pair
of sets of orientations $\bar{\Phi}_{i}^{n}(\textit{InD})$,
$\bar{\Phi}_{i}^{n}(G)$ for example. We use $\delta$, invariance score, as the
similarity metric, which is defined (based on previous work [5]) in the
following way:

$\begin{split}\delta(\bar{\Phi}_{i}^{n}(\mathit{InD}),\bar{\Phi}_{i}^{n}(G))=1-\left|\frac{\bar{\Phi}_{i}^{n}(G)-\bar{\Phi}_{i}^{n}(\mathit{InD})}{\bar{\Phi}_{i}^{n}(G)+\bar{\Phi}_{i}^{n}(\mathit{InD})}\right|.\end{split}$
(6)

We note that under some conditions, $\delta$ reports a high, yet trivial,
invariance. Namely, if the response of a neuron is low or zero for both
elements of the pair, the denominator approaches zero and the invariance
becomes large. However in this case the neuron is not responding to anything —
any activity is most likely noise. We therefore calculate a threshold of
activity for neural response invariances to be considered to contribute to the
generalization capability of the network. Otherwise, these invariances are not
integrated into the overall network invariance metric. The threshold, $\tau$,
is the 95th percentile of activity for all neurons across all images. We
employ $\tau$ with an indicator function as follows:

$\begin{split}\@ADDCLASS{ltx_eqn_lefteqn}$\displaystyle\mathbf{1}(\bar{\Phi}_{i}^{n}(\mathit{InD}),\bar{\Phi}_{i}^{n}(G))$\mbox{}\hfil\\\
&:=\begin{cases}1&\text{if
}\bar{\Phi}_{i}^{n}(\mathit{InD})\geq\tau\wedge\bar{\Phi}_{i}^{n}(G)\geq\tau\\\
0&\text{otherwise}\end{cases}.\end{split}$

Finally, we can compute the overall network generalizable and non-
generalizable invariance scores. To do so, we compute a triple average: an
average activation over the set of orientations (Eq. 5) and averaged over the
invariance of all neurons and object instances. We say that the generalizable
invariance score is the invariance between the in-distribution orientations
and the generalizable orientations determined as follows:

$\frac{1}{L}\sum_{n\in N}\sum_{i\in
I}\mathbf{1}(\bar{\Phi}^{n}_{i}(\mathit{InD}),\bar{\Phi}^{n}_{i}(G))\cdot\delta(\bar{\Phi}^{n}_{i}(\mathit{InD}),\bar{\Phi}^{n}_{i}(G)),$
(7)

where $L$ is the quantity of activity pairs above the threshold $\tau$,
_i.e.,_

$L=\sum_{n\in N}\sum_{i\in
I}\mathbf{1}(\bar{\Phi}^{n}_{i}(\mathit{InD}),\bar{\Phi}^{n}_{i}(G)).$ (8)

The definition of the network’s non-generalizable invariance score is the
same, though $\neg G$ replaces $G$.

###  4.4 Experimental Controls

a b c
---
Figure 5: Object Datasets. In our experiments we used three object categories:
(a) Airplanes, (b) Cars, and (c) Shepard&Metzler objects. The first two were
curated from ShapeNet [22] and we procedurally generated the last one. There
are $50$ instances per object category (_e.g.,_ ‘Concorde’ or ‘Spitfire’ for
the Airplanes). Images were rendered from the 3D models under fixed lighting
conditions, and the models were centered and fully contained within the image
frame. For fully-seen instances (see Fig. 1), orientations were uniformly
sampled at random using Euler angles in the range of
$-\pi\leq\alpha<\pi,-\frac{\pi}{2}\leq\beta<\frac{\pi}{2},-\pi\leq\gamma<\pi$.
For partially-seen instances, orientations were uniformly sampled from a
subset of these ranges.

#### Proportion of fully-seen instances.

We vary diversity in terms of the number of fully-seen instances $N$ between
10 (20% of the total number of instances) and 40 (80%). The remaining
instances are partially-seen. For a fair evaluation of the effect of data
diversity, the amount of training examples is kept constant as we vary the
data diversity.

#### Object Categories.

We used three categories of objects: Airplanes, Cars and Shepard&Metzler
objects. For the airplanes and cars we curated 50 high quality object
instances of each category from the ShapeNet[22] database. Both airplanes and
cars have clear axes of symmetry, which allow for intuition of how networks
generalize to OoD orientations. We therefore also experimented with highly
asymmetric objects similar to those tested for 3D mental rotations in [15]
(which we denote as Shepard&Metzler objects; Fig. 5).

#### DNN Architectures.

We used ResNet18 [14], DenseNet [23] and CORnet [24] in our experiments. The
first two were chosen as they are representative feed-forward DNNs. The
architecture of CORnet is brain-inspired and includes recurrence at higher
layers in addition to convolutions in lower layers.

#### Repetition.

We re-run each experiment five times, each time randomly sampling the specific
instances which comprise the fully-seen and partially seen sets.

#### Hyperparameters for training.

We trained the three deep convolutional neural networks using the Adam
Optimizer [25] with following learning rates and batch seizes, respectively:

  * •

ResNet18

    * –

Learning rate: 0.001

    * –

Batch size: 230

  * •

DenseNet

    * –

Learning rate: 0.001

    * –

Batch size: 64

  * •

CORnet

    * –

Learning rate: 0.0001

    * –

Batch size: 128

Batch sizes were chosen to be as large as possible while still fitting the
model, the batch of images and forward-pass computations in memory. Learning
rates were chosen from $10^{x},x\in\\{-1,-2,-3,-4,-5\\}$ to be as large as
possible while ensuring that OoD generalization remained stable. Each network
was trained for 10 epochs. After this point in-distribution performance was
stabilized at 100% and OoD performance reached an asymptote.

#### Dataset Size.

Each dataset is 200k images, 4k image for each of the 50 object instances. A
training epoch iterates through every image in the dataset once.

#### Hardware details.

Experiments were run with one CPU, 25GB of memory and on several generations
of Nvidia GPUs with a minimum of 11GB of memory.

## Author Contributions

AC, XB, DH, SM designed research; AC performed experiments with contributions
of XB; AC, XB, DH and TS analyzed data with contributions of SM; AC, XB, DH
and TS wrote the paper with contributions of SM and PS; HP and PS supervised
the research with contributions of XB, DH and TS.

## Acknowledgments

We are grateful to Tomaso Poggio and Shimon Ullman for their insightful advice
and warm encouragement. This work was supported by Fujitsu Limited (Contract
No. 40008819 and 40009105), NSF grant IIS-1901030, and by the Center for
Brains, Minds and Machines (funded by NSF STC award CCF-1231216). PS and XB
are supported by the R01EY020517 grant from the National Eye Institute (NIH).

This research is based upon work supported in part by the Office of the
Director of National Intelligence (ODNI), Intelligence Advanced Research
Projects Activity (IARPA), via [2022-21102100009]. The views and conclusions
contained herein are those of the authors and should not be interpreted as
necessarily representing the official policies, either expressed or implied,
of ODNI, IARPA, or the U.S. Government. The U.S. Government is authorized to
reproduce and distribute reprints for governmental purposes notwithstanding
any copyright annotation therein.

## Conflicts of Interests Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest. Fujitsu Limited funded this study (Contract No. 40008819
and 40009105) and also participated in the study through AC, XB and TS. All
authors declare no other competing interests.

## Data and Code Availability Statement

The raw data and code supporting the conclusions of this article are publicly
accessible via https://github.com/avicooper1/OOD_Orientation_Generalization.

## References

  * [1] Pawan Sinha and Tomaso Poggio.  Role of learning in three-dimensional form perception.  Nature, 384(6608):460–463, 1996.
  * [2] Shimon Ullman.  High-level vision: Object recognition and visual cognition.  MIT Press, 1996.
  * [3] Andrei Barbu, David Mayo, Julian Alverio, William Luo, Christopher Wang, Dan Gutfreund, Josh Tenenbaum, and Boris Katz.  ObjectNet: A large-scale bias-controlled dataset for pushing the limits of object recognition models.  In Advances in Neural Information Processing Systems, pages 9448–9458, 2019.
  * [4] Michael A Alcorn, Qi Li, Zhitao Gong, Chengfei Wang, Long Mai, Wei-Shinn Ku, and Anh Nguyen.  Strike (with) a pose: Neural networks are easily fooled by strange poses of familiar objects.  In Proc of the IEEE/CVF Conf on Computer Vision and Pattern Recognition, pages 4845–4854, 2019.
  * [5] Spandan Madan, Timothy Henry, Jamell Dozier, Helen Ho, Nishchal Bhandari, Tomotake Sasaki, Frédo Durand, Hanspeter Pfister, and Xavier Boix.  When and how convolutional neural networks generalize to out-of-distribution category–viewpoint combinations.  Nature Machine Intelligence, 4(2):146 – 153, 2022.
  * [6] Demis Hassabis, Dharshan Kumaran, Christopher Summerfield, and Matthew Botvinick.  Neuroscience-inspired artificial intelligence.  Neuron, 95(2):245–258, 2017.
  * [7] Shimon Ullman.  Using neuroscience to develop artificial intelligence.  Science, 363(6428):692–693, 2019.
  * [8] Nikos K Logothetis and Jon Pauls.  Psychophysical and physiological evidence for viewer-centered object representations in the primate.  Cerebral Cortex, 5(3):270–288, 1995.
  * [9] Robert Desimone, Thomas D Albright, Charles G Gross, and Charles Bruce.  Stimulus-selective properties of inferior temporal neurons in the macaque.  Journal of Neuroscience, 4(8):2051–2062, 1984.
  * [10] Eucaly Kobatake and Keiji Tanaka.  Neuronal selectivities to complex object features in the ventral visual pathway of the macaque cerebral cortex.  Journal of neurophysiology, 71(3):856–867, 1994.
  * [11] Nikos K Logothetis and David L Sheinberg.  Visual object recognition.  Annual review of neuroscience, 19(1):577–621, 1996.
  * [12] Tomaso Poggio and Fabio Anselmi.  Visual cortex and deep networks: learning invariant representations.  MIT Press, 2016.
  * [13] Syed Suleman Abbas Zaidi, Xavier Boix, Neeraj Prasad, Sharon Gilad-Gutnick, Shlomit Ben-Ami, and Pawan Sinha.  Is robustness to transformations driven by invariant neural representations?  arXiv preprint, arXiv:2007.00112, 2020.
  * [14] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proc of the IEEE/CVF Conf on Computer Vision and Pattern Recognition, pages 770–778, 2016.  We used the following implementation in our experiments: https://pytorch.org/vision/stable/models.html#torchvision.models.resnet18.
  * [15] Roger N. Shepard and Jacqueline Metzler.  Mental rotation of three-dimensional objects.  Science, 171(3972):701–703, 1971.
  * [16] Wang Angtian, Adam Kortylewski, and Alan Yuille.  Nemo: Neural mesh models of contrastive features for robust 3d pose estimation.  In Proc of the Int Conf on Learning Representations, 2021.
  * [17] Taco S Cohen, Mario Geiger, Jonas Köhler, and Max Welling.  Spherical CNNs.  In Proc of the Int Conf on Learning Representations, 2018.
  * [18] Holly A Ruff.  Effect of object movement on infants’ decision of object structure.  Developmental Psychology, 18(3):462, 1982.
  * [19] Scott P Johnson and Richard N Aslin.  Perception of object unity in young infants: The roles of motion, depth, and orientation.  Cognitive Development, 11(2):161–180, 1996.
  * [20] Daniel LK Yamins and James J DiCarlo.  Using goal-driven deep learning models to understand sensory cortex.  Nature neuroscience, 2016.
  * [21] Herbert Goldstein, Charles Poole, and John Safko.  Classical mechanics.  Addison-Wesley, 3rd edition, 2002.
  * [22] Angel X. Chang, Thomas Funkhouser, Leonidas Guibas, Pat Hanrahan, Qixing Huang, Zimo Li, Silvio Savarese, Manolis Savva, Shuran Song, Hao Su, Jianxiong Xiao, Li Yi, and Fisher Yu.  ShapeNet: An Information-Rich 3D Model Repository.  Technical Report arXiv:1512.03012 [cs.GR], Stanford University — Princeton University — Toyota Technological Institute at Chicago, 2015.
  * [23] Gao Huang, Zhuang Liu, Laurens Van Der Maaten, and Kilian Q Weinberger.  Densely connected convolutional networks.  In Proc of the IEEE/CVF Conf on Computer Vision and Pattern Recognition, pages 4700–4708, 2017.  We used the following implementation in our experiments: https://pytorch.org/vision/stable/models.html#torchvision.models.densenet121.
  * [24] Jonas Kubilius, Martin Schrimpf, Aran Nayebi, Daniel Bear, Daniel L. K. Yamins, and James J. DiCarlo.  CORnet: Modeling the neural mechanisms of core object recognition.  bioRxiv preprint, 408385, 2018.  We used the following implementation in our experiments: https://github.com/dicarlolab/CORnet.
  * [25] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint, arXiv:1412.6980, 2017.
  * [26] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, et al.  ImageNet Large Scale Visual Recognition Challenge.  International Journal of Computer Vision, 115(3):211–252, 2015\.

## Supplementary Figures

a

|  |
---|---|---

b

|  |
---|---|---

Figure S.1: Accuracy heatmaps: alternative ’seed’ orientations. (a) ’Seed’
(in-distribution) orientations include $-0.25\leq\alpha\leq
0.25,-0.25\leq\gamma\leq 0.25,-1/2\pi\leq\beta<1/2\pi$. (b) ’Seed’
orientations include $-0.1\leq\beta\leq 0.1,-0.25\leq\gamma\leq
0.25,-1.8\pi\leq\alpha<-1.3/\pi$.

a

|  |
---|---|---
20 Fully Seen |  30 Fully Seen |  40 Fully Seen

b

|  |
---|---|---
20 Fully Seen |  30 Fully Seen |  40 Fully Seen

Figure S.2: Accuracy heatmaps: effect of data diversity - alternative object
categories. Increasing number of fully-seen instances, with different object
classes. (a) Shepard-Metzler Objects. (b) Cars.

a

|  |
---|---|---

b

|  |
---|---|---

Figure S.3: Accuracy heatmaps: alterative training conditions - pretraining
and augmentation. (a) ResNet-18 pretrained on ImageNet [26], finetuned on our
learning paradigm with airplanes. Network behavior isn’t meaingfully altered.
(b) All data (both from fully-seen and partially-seen instances) were
augmented with random 2D image rotations. This effectively expands the in-
distribution set to include all generalizable orientations. This results in
generalizable orientations with high accuracy.

a

|  |
---|---|---

b

|  |
---|---|---

Figure S.4: Accuracy heatmaps: alternative backbone architectures. Network’s
backbone used (in place of ResNet-18): (a) DenseNet. (b) CORnet.

Figure S.5: Modeling generalization patterns for OoD orientations, continued.
The same analysis as Figs. 3b is applied to the controls introduced in Figs.
S.3, S.4.

OoD Accuracy, ResNet, Base Experiments
OoD Accuracy, Various Controls

Figure S.6: OoD accuracy, split between generalizable and non-generalizable
orientations In Fig. 3a we report the average accuracy across all OoD
orientations. As we note, however, accuracy behavior is differentiated between
generalizable and non-generalizable orientations. Here we report the average
accuracy for these two orientation groups. Gray horizontal lines indicate
chance performance of 2% and 10% (the latter relevant in the case where fully-
seen and partially-seen instances are of two different classes.) Generalizable
accuracy is always greater than non-generalizable accuracy. The former is
always well above chance, while the latter is below or at chance level. (a)
The generalizable and non-generalizable average accuracy for the same set of
experiments presented in Fig. 3a. (b) The average accuracies for several other
conditions. These other conditions are explained in Figs. S.3, S.4.

a | b
---|---
|
Figure S.7: Invariance and Dissemination: controls. The same analysis as Figs.
4b, c is applied to the controls introduced in Figs. S.3, S.4.