provides very competitive results, showing the utility of our approach in
dealing with unobserved sources of heterogeneity in the treatment response.

The performance of CDVAE is very close to that of the best model (CRN), while
the latter takes as an additional input the cluster labels. This shows that
despite the lack of any information related to the subgroups of individuals,
CDVAE manages to provide a very competitive prediction of the counterfactual
response, confirming further the reliability of the constructed substitutes.

##  10 Conclusion

Our work introduces CDVAE, a novel model for estimating Individual Treatment
Effects (ITE) in longitudinal data. CDVAE combines Dynamic Variational
Autoencoder (DVAE) techniques to address the challenges posed by unobserved
adjustment variables and time-varying effects. We establish the theoretical
validity of CDVAE, ensuring the learned adjustment variables are a reliable
substitute for unobserved factors. Through extensive experiments, CDVAE
outperforms state-of-the-art models in estimating ITE, expanding the
possibilities of causal inference in longitudinal settings. Our work focused
on a binary treatment; thus, a natural extension considers continuous
treatments, possibly multiple. Recent works attempt to derive generalization
bounds over the conditional treatment effect in such a setting but only in the
static case (Wang et al., 2022; Bellot et al., 2022). The same lack of
theoretical guarantees applies to the effect of a sequence of interventions in
a dynamic regime. Another perspective is to go beyond estimating ITEs for
panel data with fixed time length and perform counterfactual forecasting of
ITEs over multiple time steps ahead, given the same hypothesis on
heterogeneity. Finally, we can be interested in inferring valid substitutes
for time-varying risk factors, i.e., within-individual sources of variation as
explained in Subsection 3.1. Inspiration may be possible from works on
inferring time-varying confounders when sequential ignorability is violated
(Bica et al., 2020b; Hatt and Feuerriegel, 2021).

##  11 Appendix A. Algorithmic Details

In this appendix, we provide algorithmic details for CDVAE. First, we compute
an approximation of the Integral Probability Metric (IPM) term used in (32).
The IPM term aims to reduce the covariate imbalance between the treatment and
control groups by quantifying the dissimilarity between their distributions.
We utilize the Sinkhorn-Knopp algorithm (Sinkhorn, 1967) and the Wasserstein
distance computation algorithm (algorithm 3 in Cuturi and Doucet (2014)) to
calculate the IPM. The Wasserstein distance computation (algorithm 2)
calculates the pairwise distances between data points and constructs a kernel
matrix using a regularization parameter. It then computes the row and column
marginals based on the weights assigned to each data point. The algorithm
utilizes the Sinkhorn-Knopp algorithm (algorithm 1), included within the
Wasserstein distance computation algorithm, to compute the optimal transport
matrix. Finally, the Wasserstein distance is obtained by summing the products
of the optimal transport matrix and the pairwise distances.

Algorithm 1 Sinkhorn-Knopp Algorithm

1:Kernel matrix $K\in\mathbb{R}^{n_{t}\times n_{c}}$.

2:Row marginal vector $a\in\mathbb{R}^{n_{t}}$, Column marginal vector
$b\in\mathbb{R}^{n_{c}}.$

3:Optimal transport matrix $T\in\mathbb{R}^{n_{t}\times n_{c}}.$

4:Initialize transport matrix $T^{(0)}$ with all entries set to $1$.

5:Set iteration counter $k=0$.

6:while not converged do

7: Update row scaling vector $u\in\mathbb{R}^{n_{t}}$.

8: $u_{i}\leftarrow\frac{a_{i}}{\sum_{j=1}^{n_{c}}K_{ij}T^{(k)}_{ij}}$.

9: Update column scaling vector $v\in\mathbb{R}^{n_{c}}$.

10: $v_{j}\leftarrow\frac{b_{j}}{\sum_{i=1}^{n_{t}}K_{ij}T^{(k)}_{ij}}$.

11: Update transport matrix $T^{(k+1)}$.

12:
$T^{(k+1)}_{ij}\leftarrow\frac{u_{i}K_{ij}v_{j}}{\sum_{i^{\prime}=1}^{n_{t}}\sum_{j^{\prime}=1}^{n_{c}}u_{i^{\prime}}K_{i^{\prime}j^{\prime}}v_{j^{\prime}}}$.

13: Increment iteration counter $k\leftarrow k+1$.

14: if convergence criterion met then

15: Return $T^{(k+1)}$.

16: end if

17:end while

Algorithm 2 Weighted Wasserstein Distance Computation

1:Batch
$\mathcal{B}=\\{\\{\omega_{it},y_{it},\mathbf{x}_{it}\\}_{t=1}^{T},i=1,\dots,|\mathcal{B}|\\}$.

2:Representation learner $\Phi$.

3:Weights vectors $\alpha_{\Phi}(\mathbf{h}_{t},\omega)$.

4:Regularization parameter $\lambda$.

5:for $t\in\\{1,2,\dots,T\\}$ do

6: Compute $n_{t}^{(t)}=\sum_{i\in\mathcal{B}}W_{it}$,
$n_{c}^{(t)}=\sum_{i\in\mathcal{B}}1-W_{it}$.

7: Compute pairwise distances matrix $M\in\mathbb{R}^{n_{t}^{(t)}\times
n_{c}^{(t)}}$ {$M_{ij}^{(t)}=\|\mathbf{H}_{it}-\mathbf{H}_{jt}\|_{L_{2}}$,
$\forall i,j\in\mathcal{B}$; $W_{it}=1$ and $W_{jt}=0$.}

8: Initialize kernel matrix $K\in\mathbb{R}^{n_{t}^{(t)}\times n_{c}^{(t)}}$
such as $K_{ij}^{(t)}=e^{-\lambda M_{ij}^{(t)}}$.

9: Compute row marginal vector $a^{(t)}\in\mathbb{R}^{n_{t}}$ such as
$a_{i}^{(t)}=\frac{\alpha_{\Phi}(\mathbf{h}_{it},1)}{\sum_{k=1,W_{kt}=1}\alpha_{\Phi}(\mathbf{h}_{it},1)}$

10: Compute column marginal vector $b^{(t)}\in\mathbb{R}^{n_{c}}$ such as
$b_{j}^{(t)}=\frac{\alpha_{\Phi}(\mathbf{h}_{it},0)}{\sum_{k=1,W_{kt}=0}\alpha_{\Phi}(\mathbf{h}_{it},0)}$

11: Compute optimal transport matrix $T^{(t)}\in\mathbb{R}^{n_{t}^{(t)}\times
n_{c}^{(t)}}$: $T\leftarrow\text{Sinkhorn-Knopp}(K^{(t)},a^{(t)},b^{(t)})$.

12: Compute Wasserstein distance
$D_{t}=\sum_{i=1}^{n_{t}}\sum_{j=1}^{n_{c}}T_{ij}^{(t)}M_{ij}^{(t)}$.

13:end for

14:Return $\sum_{t=1}^{T}D_{t}$.

Algorithm 3 Pseudo-code for training $\beta_{cyc}$-CDVAE

1:Training Data
$\mathcal{D}_{T}=\\{\\{w_{it},y_{it},\mathbf{x}_{it}\\}_{t=1}^{T},i=1,\dots,n\\}$.

2:CDVAE parameters ($\phi$ , $\theta_{y}$, $\theta_{\omega}$, parameters of
shared representation $\Phi$), Optimizer parameters.

3:for $p\in\\{1,\dots,\mathrm{epoch_{max}}\\}$ do

4: for
$\mathcal{B}=\\{\\{w_{it},y_{it},\mathbf{x}_{it}\\}_{t=1}^{T},i=1,\dots,|\mathcal{B}|\\}$
do

5: Compute approximate posterior $q_{\phi}(z|y_{\leq T},\mathbf{x}_{\leq
T},\omega_{\leq T})$.

6: Sample latent variables $z$ from $q_{\phi}(z|y_{\leq T},\mathbf{x}_{\leq
T},\omega_{\leq T})$.

7: Compute reconstruction $\hat{y}_{t}$ using the decoder network.

8: Compute weighted reconstruction loss for $\hat{y}_{t}$ and $y_{t}$.

9: Compute KL divergence loss between $q_{\phi}(z|y_{\leq T},\mathbf{x}_{\leq
T},\omega_{\leq T})$ and prior $p(z)$.

10: Compute IPM term using algorithm 2.

11: Compute moment matching penalty.

12: Compute binary cross-entropy loss for the propensity network.

13: Update $\beta_{cyc}$ according to the current iteration number.

14: Update $\beta_{cyc}$-CDVAE parameters using gradients from total loss
(32).

15: end for

16:end for

17:Return Trained CDVAE model.

##  12 Appendix B. Experimental Details

###  12.1 Training Strategy

For each simulation of a synthetic data set, we generate 10000 units and
divide them randomly into 70% for training data, 15% for validation, and 15%
for out-of-sample testing. We use Pytorch (Paszke et al., 2019) and Pytorch
Lightening (Falcon and The PyTorch Lightning team, 2019) to implement CDVAE,
RMSMs, and CRN, while we use the EconML (Battocchi et al., 2019) package for
the CausalForestDML. For a single model training, we fix the batch size to 128
and the maximal number of epochs to 100. We use stop loss to end training once
the validation loss increases or stops improving after four consecutive
epochs. We fix the weight-decay to $10^{-4}$ and use the Stochastic Average
Weighting (SWA) strategy (Izmailov et al., 2018) during training to improve
generalization performance. The learning rate of SWA is set to $10^{-2}$ and
the number of annealing epochs to 3. We clip gradients’ global norm to less
than $0.5$ to avoid exploding gradients during optimization. The cycling
strategy was fixed during all the experiments: the maximal number of cycles is
$M=6$. At each cycle, the parameter $\beta$ increases linearly from 0 to 1
half of the time, i.e., $R=0.5$, and then stagnates at 1 for the second half
of the cycle (cf. Figure 12).

Figure 12: Cyclical annealing of the hyperparameter $\beta$ during training of
$N_{iter}=12000$ maximum iterations. The parameter $\beta$ evolves following
$M=6$ cycles. Each cycle starts with a monotone increase (linear) of the
penalty until it reaches the value $\beta=1$ corresponding to the standard,
non-penalized ELBO.

For the hyperparameters selection, we fine-tune our models using the Bayesian
optimization technique, namely the TPE algorithm (Bergstra et al., 2011) using
the Optuna library (Akiba et al., 2019). We use the weighted reconstruction
error as the selection criterion (objective function for the Bayesian
optimization). We do not use a metric related to the quality of estimating ITE
as a selection criterion; a choice cannot be used in a real scenario because
ITEs are usually unavailable. It is still rather an open problem how to design
a criterion for causal cross-validation or hyperparameters tuning for a causal
model (Dorie et al., 2019). Many proxies are used in the literature, including
loss over the factual outcomes (Lim, 2018; Bica et al., 2020a; Hassanpour and
Greiner, 2019a; Bica et al., 2020b), one nearest-neighbor imputation (Shalit
et al., 2017; Johansson et al., 2022), for the counterfactual outcome,
influence functions (Alaa and Van Der Schaar, 2019), rank-preserving causal
cross-validation (Schuler et al., 2018), Robinson residual decomposition (Nie
and Wager, 2021; Lu et al., 2020). In this work, and since we intend to search
for the best regularization parameters ($\lambda_{IPM}$, $\lambda_{MM}$,
$\lambda_{W}$) among other parameters related to CDVAE architecture and
optimization, we chose not to consider the total loss as this may bias the
choice of ( $\lambda_{IPM}$, $\lambda_{MM}$, $\lambda_{W}$) toward small
values. We, therefore, use the weighted reconstruction error as an objective
function for fine-tuning.

Similarly, we fine-tune CRN and RMSM using the loss over the factual response
as a criterion. We implement fine-tuning of the first stage modeling in
CausalForestDML bet setting the criterion for the treatment model (a Gradient
Boosting Classifier) to be the accuracy and for the outcome model (a Gradient
Boosting Regressor) to be the R-score. The second stage of modeling in
CausalForestDML was tuned using a GridSearch, whose parameters and search
space are internally set in EconML.

We report in Table 5 the search space of hyperparameters for CDVAE and in
Table6 the best-found hyperparameters values. Similarly, we report the best
hyperparameters values found for CRN (Tables 7 and 8), RMSM (Tables 10, 9, 11,
and 12), and first stage estimation in CausalForestDML (Tables 13 and 14).

Hyperparameter | Search Space
---|---
Hidden size (LSTM) | $[16,18,20,\ldots,32,34,36]$
Dimension of $\phi(\mathbf{h}_{t})$ | $[16,18,20,\ldots,32,34,36]$
Dimension of $\mathbf{z}$ | $[16,18,20,\ldots,32,34,36]$
Learning rate | Log-uniform distribution between $1\times 10^{-5}$ and $1\times 10^{-2}$
LSTM layers | Integer values from 1 to 3
LSTM dropout probability | Uniform distribution between $0.0$ and $0.5$
$\lambda_{W}$ | Uniform distribution between $0.0$ and $5.0$
$\lambda_{ipm}$ | Uniform distribution between $0.0$ and $5.0$
$\lambda_{mm}$ | Uniform distribution between $0.0$ and $5.0$
Table 5: Search Space for CDVAE Hyperparameters Hyperparameter | Value
---|---
Hidden size (LSTM) | 32
Dimension of $\phi(\mathbf{h}_{t})$ | 36
Dimension of $\mathbf{z}$ | 22
Learning rate | 0.0003
LSTM layers | 1
LSTM dropout probability | -
$\lambda_{W}$ | 0.65
$\lambda_{ipm}$ | 0.45
$\lambda_{mm}$ | 3.75
Table 6: Best Hyperparameters for CDVAE Hyperparameter | Description
---|---
hidden size LSTM | 32
latent_dim | 36
LSTM layers | 2
LSTM dropout probability | 0.15
learning_rate | 0.0002
Table 7: Hyperparameters of the CRN encoder model Hyperparameter | Description
---|---
Decoder latent dimension | 16
LSTM dropout probability | -
Learning rate | 0.0035
LSTM layers | 1
Table 8: Hyperparameters of the CRN decoder model Hyperparameter | Description
---|---
hidden size LSTM | 16
latent_dim | 36
LSTM dropout probability | -
Learning rate | 0.009
LSTM layers | 1
Table 9: Hyperparameters of the RMSM propensity history network Hyperparameter | Description
---|---
hidden size LSTM | 32
latent_dim | 20
LSTM dropout probability | -
Learning rate | 0.0035
LSTM layers | 1
Table 10: Hyperparameters of the RMSM propensity network Hyperparameter | Description
---|---
hidden size LSTM | 32
LSTM dropout probability | 0.4
Learning rate | 0.0002
LSTM layers | 2
Table 11: Hyperparameters of the RMSM encoder Hyperparameter | Description
---|---
Decoder latent dimension | 24
LSTM dropout probability | -
Learning rate | 0.0003
LSTM layers | 1
Table 12: Hyperparameters of the RMSM decoder Hyperparameter | Value
---|---
Learning rate | 0.016
Max depth | 9
Min samples leaf | 7
min samples split | 7
Number of estimators | 750
Table 13: Hyperparameters of the treatment model (CausalForestDML) Hyperparameter | Value
---|---
Learning rate | 0.015
Max depth | 4
Min samples leaf | 8
min samples split | 5
Number of estimators | 350
Table 14: Hyperparameters of the outcome model (CausalForestDML)

###  12.2 The Extended Neural Architecture of CDVAE

The extended neural architecture of CDVAE comprises multiple components which
we did not explicit in Section 8. We first begin by detailing neural network
functions related to the generative model. The Table outlines the architecture
for the Representation Learner $\Phi$ which encodes the context history, Table
17 presents the identical architecture for both $f_{\theta_{y}^{1}}$ and
$f_{\theta_{y}^{0}}$ responsible for generating the two potential outcomes.
Meanwhile, Table 17 illustrates the design of propensity network
$e_{\theta_{\omega}}(.)$ built on the top of the shared representation.
Lastly, Tables 19 and 19 depict the architecture used to learn both the mean
and covariance matrix for the approximate posterior assumed to be Gaussian.

Inputs: $\\{y_{t},\mathbf{x}_{t},\omega_{t}\\}_{1\leq t\leq T}$
---
Concat: $[y_{t},\mathbf{x}_{t},\omega_{t}]_{1\leq t\leq T}$
LSTM layer
Linear Layer
LeakyReLU
Linear Layer
Outputs: $\\{\Phi(\textbf{h}_{t})\\}_{1\leq t\leq T}$
Table 15: Architecture for representation learner $\phi$ of CDVAE

Inputs: $\\{\Phi(\mathbf{h}_{t})\\}_{1\leq t\leq T}$, $\mathbf{z}$
---
Concatenate: $[\Phi(\mathbf{h}_{t}),\mathbf{z}]_{1\leq t\leq T}$
Linear Layer
LeakyReLU
Linear Layer
Output: $\\{\hat{Y}_{t+1}(\omega)\\}_{1\leq t\leq T-1}$
Table 16: $f_{\theta_{y}^{\omega}}$ from (34)

Inputs: $\\{\Phi(\mathbf{h}_{t})\\}_{1\leq t\leq T}$
---
Linear Layer
LeakyReLU
Linear Layer
Sigmoid
Output: $\\{\hat{W}_{t+1}(\omega)\\}_{1\leq t\leq T-1}$
Table 17: $e_{\theta_{\omega}}(.)$ from (35)

Inputs: $\mathbf{g}_{T}$
---
Linear Layer
LeakyReLU
Linear Layer
Output: $\mu_{\phi_{2}}\left(\mathbf{g}_{T}\right)$
Table 18: Mean of $q_{\phi}\left(\mathbf{z}\mid y_{\leq T},\mathbf{x}_{\leq
T},\omega_{\leq T}\right)$

Inputs: $\mathbf{g}_{T}$
---
Linear Layer
LeakyReLU
Linear Layer
Output: Diagonal $\Sigma_{\phi_{3}}\left(\mathbf{g}_{T}\right)$
Table 19: Covariance of $q_{\phi}\left(\mathbf{z}\mid y_{\leq
T},\mathbf{x}_{\leq T},\omega_{\leq T}\right)$

Acknowledgments and Disclosure of Funding

This work was funded by the French Ministry of Higher Education and Research
and Saint-Gobain under the Grant CIFRE n°2021/0476.

## References

  * Akiba et al. (2019) T. Akiba, S. Sano, T. Yanase, T. Ohta, and M. Koyama.  Optuna: A next-generation hyperparameter optimization framework.  In _Proceedings of the 25th ACM SIGKDD international conference on knowledge discovery & data mining_, pages 2623–2631, 2019.
  * Alaa and Van Der Schaar (2019) A. Alaa and M. Van Der Schaar.  Validating causal inference models via influence functions.  In _International Conference on Machine Learning_ , pages 191–201. PMLR, 2019.
  * Assaad et al. (2021) S. Assaad, S. Zeng, C. Tao, S. Datta, N. Mehta, R. Henao, F. Li, and L. Carin.  Counterfactual representation learning with balancing weights.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1972–1980. PMLR, 2021.
  * Atan et al. (2018) O. Atan, J. Jordon, and M. Van der Schaar.  Deep-treat: Learning optimal personalized treatments from observational data using neural networks.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 32, 2018.
  * Athey et al. (2019) S. Athey, J. Tibshirani, and S. Wager.  Generalized random forests.  _The Annals of Statistics_ , 2019.
  * Austin (2011) P. C. Austin.  An introduction to propensity score methods for reducing the effects of confounding in observational studies.  _Multivariate behavioral research_ , 46(3):399–424, 2011.
  * Battocchi et al. (2019) K. Battocchi, E. Dillon, M. Hei, G. Lewis, P. Oka, M. Oprescu, and V. Syrgkanis.  EconML: A Python Package for ML-Based Heterogeneous Treatment Effects Estimation.  https://github.com/microsoft/EconML, 2019.  Version 0.x.
  * Bayer and Osendorfer (2014) J. Bayer and C. Osendorfer.  Learning stochastic recurrent networks.  _ArXiv_ , abs/1411.7610, 2014.
  * Bellot et al. (2022) A. Bellot, A. Dhir, and G. Prando.  Generalization bounds and algorithms for estimating conditional average treatment effect of dosage.  _arXiv preprint arXiv:2205.14692_ , 2022.
  * Bergstra et al. (2011) J. Bergstra, R. Bardenet, Y. Bengio, and B. Kégl.  Algorithms for hyper-parameter optimization.  _Advances in neural information processing systems_ , 24, 2011.
  * Bica et al. (2020a) I. Bica, A. M. Alaa, J. Jordon, and M. van der Schaar.  Estimating counterfactual treatment outcomes over time through adversarially balanced representations.  _arXiv preprint arXiv:2002.04083_ , 2020a.
  * Bica et al. (2020b) I. Bica, A. M. Alaa, and M. van der Schaar.  Time series deconfounder: Estimating treatment effects over time in the presence of hidden confounders.  In _International Conference on Machine Learning_ , 2020b.
  * Bray et al. (2018) N. Bray, T. Doherty, and M. Montero-Odasso.  The effect of high dose vitamin d3 on physical performance in frail older adults. a feasibility study.  _The Journal of Frailty & Aging_, 7:155–161, 2018.
  * Breiman (2001) L. Breiman.  Random forests.  _Machine Learning_ , 45:5–32, 2001.
  * Cai et al. (2022) Y. Cai, A. A. Wanigatunga, C. M. Mitchell, J. K. Urbanek, E. R. Miller, S. P. Juraschek, E. D. Michos, R. R. Kalyani, D. L. Roth, L. J. Appel, et al.  The effects of vitamin d supplementation on frailty in older adults at risk for falls.  _BMC geriatrics_ , 22(1):1–9, 2022.
  * Cheng et al. (2021) L. Cheng, R. Guo, and H. Liu.  Causal mediation analysis with hidden confounders.  _Proceedings of the Fifteenth ACM International Conference on Web Search and Data Mining_ , null:null, 2021.  doi: 10.1145/3488560.3498407.  URL https://www.semanticscholar.org/paper/81209a11b05b4bd3ca254d6b0519a274d1f79723.
  * Chernozhukov et al. (2016) V. Chernozhukov, D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins.  Double/debiased machine learning for treatment and causal parameters.  _arXiv preprint arXiv:1608.00060_ , 2016.
  * Chernozhukov et al. (2017) V. Chernozhukov, M. Goldman, V. Semenova, and M. Taddy.  Orthogonal machine learning for demand estimation: High dimensional causal inference in dynamic panels.  _arXiv_ , pages arXiv–1712, 2017.
  * Conover (1999) W. J. Conover.  _Practical nonparametric statistics_ , volume 350.  john wiley & sons, 1999.
  * Crump et al. (2009) R. K. Crump, V. J. Hotz, G. W. Imbens, and O. A. Mitnik.  Dealing with limited overlap in estimation of average treatment effects.  _Biometrika_ , 96(1):187–199, 2009.
  * Cuturi and Doucet (2014) M. Cuturi and A. Doucet.  Fast computation of wasserstein barycenters.  In _International conference on machine learning_ , pages 685–693. PMLR, 2014.
  * Demidenko (2013) E. Demidenko.  _Mixed models: theory and applications with R_.  John Wiley & Sons, 2013.
  * Dorie et al. (2019) V. Dorie, J. Hill, U. Shalit, M. Scott, and D. Cervone.  Automated versus do-it-yourself methods for causal inference: Lessons learned from a data analysis competition.  2019\.
  * Falcon and The PyTorch Lightning team (2019) W. Falcon and The PyTorch Lightning team.  PyTorch Lightning, Mar. 2019.  URL https://github.com/Lightning-AI/lightning.
  * Fang et al. (2023) J. Fang, Q. Cui, G. Zhang, C. Tang, L. Gu, L. Li, J. Gu, J. Zhou, and F. Wu.  Alleviating matching bias in marketing recommendations.  In _Proceedings of the 46th International ACM SIGIR Conference on Research and Development in Information Retrieval_ , pages 3359–3363, 2023.
  * Fernández-Loría and Provost (2022) C. Fernández-Loría and F. Provost.  Causal classification: Treatment effect estimation vs. outcome prediction.  _The Journal of Machine Learning Research_ , 23(1):2573–2607, 2022.
  * Fitzmaurice et al. (2012) G. M. Fitzmaurice, N. M. Laird, and J. H. Ware.  _Applied longitudinal analysis_ , volume 998.  John Wiley & Sons, 2012.
  * Foster and Syrgkanis (2019) D. J. Foster and V. Syrgkanis.  Orthogonal statistical learning.  _ArXiv_ , abs/1901.09036, 2019.
  * Freedman and Berk (2008) D. A. Freedman and R. A. Berk.  Weighting regressions by propensity scores.  _Evaluation review_ , 32(4):392–409, 2008.
  * Fu et al. (2019) H. Fu, C. Li, X. Liu, J. Gao, A. Celikyilmaz, and L. Carin.  Cyclical annealing schedule: A simple approach to mitigating kl vanishing.  _arXiv preprint arXiv:1903.10145_ , 2019.
  * Geng et al. (2017) C. Geng, H. Paganetti, and C. Grassberger.  Prediction of treatment response for combined chemo-and radiation therapy for non-small cell lung cancer patients using a bio-mathematical model.  _Scientific reports_ , 7(1):1–12, 2017.
  * Girin et al. (2021) L. Girin, S. Leglaive, X. Bie, J. Diard, T. Hueber, and X. Alameda-Pineda.  Dynamical variational autoencoders: A comprehensive review.  _ArXiv_ , abs/2008.12595, 2021.
  * Gretton et al. (2012) A. Gretton, K. M. Borgwardt, M. J. Rasch, B. Schölkopf, and A. Smola.  A kernel two-sample test.  _The Journal of Machine Learning Research_ , 13(1):723–773, 2012.
  * Hainmueller (2012) J. Hainmueller.  Entropy balancing for causal effects: A multivariate reweighting method to produce balanced samples in observational studies.  _Political Analysis_ , 20:25 – 46, 2012.
  * Hair Jr and Sarstedt (2021) J. F. Hair Jr and M. Sarstedt.  Data, measurement, and causal inferences in machine learning: opportunities and challenges for marketing.  _Journal of Marketing Theory and Practice_ , 29(1):65–77, 2021.
  * Hansen (1982) L. P. Hansen.  Large sample properties of generalized method of moments estimators.  _Econometrica_ , 50:1029–1054, 1982.
  * Hassanpour and Greiner (2019a) N. Hassanpour and R. Greiner.  Learning disentangled representations for counterfactual regression.  In _International Conference on Learning Representations_ , 2019a.
  * Hassanpour and Greiner (2019b) N. Hassanpour and R. Greiner.  Counterfactual regression with importance sampling weights.  In _IJCAI_ , pages 5880–5887, 2019b.
  * Hatt and Feuerriegel (2021) T. Hatt and S. Feuerriegel.  Sequential deconfounding for causal inference with unobserved confounders.  _arXiv preprint arXiv:2104.09323_ , 2021.
  * Hernán MA (2020) R. J. Hernán MA.  _Causal Inference: What If_.  Boca Raton: Chapman & Hall/CRC, 2020.
  * Hill (2011) J. L. Hill.  Bayesian nonparametric modeling for causal inference.  _Journal of Computational and Graphical Statistics_ , 20(1):217–240, 2011.
  * Hochreiter and Schmidhuber (1997) S. Hochreiter and J. Schmidhuber.  Long short-term memory.  _Neural computation_ , 9(8):1735–1780, 1997.
  * Holland (1986) P. W. Holland.  Statistics and causal inference.  _Journal of the American statistical Association_ , 81(396):945–960, 1986.
  * Imai and Ratkovic (2014) K. Imai and M. T. Ratkovic.  Covariate balancing propensity score.  _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 76, 2014.
  * Imai and Ratkovic (2015) K. Imai and M. T. Ratkovic.  Robust estimation of inverse probability weights for marginal structural models.  _Journal of the American Statistical Association_ , 110:1013 – 1023, 2015.
  * Imbens and Rubin (2015) G. W. Imbens and D. B. Rubin.  _Causal inference in statistics, social, and biomedical sciences_.  Cambridge University Press, 2015.
  * Izmailov et al. (2018) P. Izmailov, D. Podoprikhin, T. Garipov, D. P. Vetrov, and A. G. Wilson.  Averaging weights leads to wider optima and better generalization.  _ArXiv_ , abs/1803.05407, 2018.
  * Johansson et al. (2016) F. Johansson, U. Shalit, and D. Sontag.  Learning representations for counterfactual inference.  In _International conference on machine learning_ , pages 3020–3029. PMLR, 2016.
  * Johansson et al. (2018) F. D. Johansson, N. Kallus, U. Shalit, and D. Sontag.  Learning weighted representations for generalization across designs.  _arXiv preprint arXiv:1802.08598_ , 2018.
  * Johansson et al. (2019) F. D. Johansson, D. Sontag, and R. Ranganath.  Support and invertibility in domain-invariant representations.  In _The 22nd International Conference on Artificial Intelligence and Statistics_ , pages 527–536. PMLR, 2019.
  * Johansson et al. (2022) F. D. Johansson, U. Shalit, N. Kallus, and D. Sontag.  Generalization bounds and representation learning for estimation of potential outcomes and causal effects.  _The Journal of Machine Learning Research_ , 23(1):7489–7538, 2022.
  * Jung et al. (2020) Y. Jung, J. Tian, and E. Bareinboim.  Learning causal effects via weighted empirical risk minimization.  _Advances in neural information processing systems_ , 33:12697–12709, 2020.
  * Kallenberg (2021) O. Kallenberg.  Foundations of modern probability.  _Probability Theory and Stochastic Modelling_ , 2021.
  * Kallus (2020) N. Kallus.  Deepmatch: Balancing deep covariate representations for causal inference using adversarial training.  In _International Conference on Machine Learning_ , pages 5067–5077. PMLR, 2020.
  * Krishnan et al. (2015) R. G. Krishnan, U. Shalit, and D. Sontag.  Deep kalman filters.  _arXiv preprint arXiv:1511.05121_ , 2015.
  * Kuroki and Pearl (2014) M. Kuroki and J. Pearl.  Measurement bias and effect restoration in causal inference.  _Biometrika_ , 101:423–437, 2014.  doi: 10.1093/BIOMET/AST066.  URL https://www.semanticscholar.org/paper/19dabe9fdbe82d74a584f5ece882f5825e4effdb.
  * LeCun et al. (2015) Y. LeCun, Y. Bengio, and G. Hinton.  Deep learning.  _nature_ , 521(7553):436–444, 2015.
  * Li et al. (2018) F. Li, K. L. Morgan, and A. M. Zaslavsky.  Balancing covariates via propensity score weighting.  _Journal of the American Statistical Association_ , 113(521):390–400, 2018.
  * Li and Greene (2013) L. Li and T. Greene.  A weighting analogue to pair matching in propensity score analysis.  _The international journal of biostatistics_ , 9(2):215–234, 2013.
  * Li et al. (2021) R. Li, S. Hu, M. Lu, Y. Utsumi, P. Chakraborty, D. M. Sow, P. Madan, J. Li, M. F. Ghalwash, Z. Shahn, and L. wei H. Lehman.  G-net: a recurrent network approach to g-computation for counterfactual prediction under a dynamic treatment regime.  In _ML4H@NeurIPS_ , 2021.
  * Li and Mandt (2018) Y. Li and S. Mandt.  Disentangled sequential autoencoder.  In _International Conference on Machine Learning_ , 2018.
  * Lim (2018) B. Lim.  Forecasting treatment responses over time using recurrent marginal structural networks.  _advances in neural information processing systems_ , 31, 2018.
  * Lopez and Gutman (2017) M. J. Lopez and R. Gutman.  Estimation of causal effects with multiple treatments: a review and new ideas.  _Statistical Science_ , pages 432–454, 2017.
  * Louizos et al. (2017) C. Louizos, U. Shalit, J. M. Mooij, D. Sontag, R. Zemel, and M. Welling.  Causal effect inference with deep latent-variable models.  _Advances in neural information processing systems_ , 30, 2017.
  * Lu et al. (2020) D. Lu, C. Tao, J. Chen, F. Li, F. Guo, and L. Carin.  Reconsidering generative objectives for counterfactual reasoning.  _Advances in Neural Information Processing Systems_ , 33:21539–21553, 2020.
  * Lunceford and Davidian (2004) J. K. Lunceford and M. Davidian.  Stratification and weighting via the propensity score in estimation of causal treatment effects: a comparative study.  _Statistics in medicine_ , 23(19):2937–2960, 2004.
  * Melnychuk et al. (2022) V. Melnychuk, D. Frauen, and S. Feuerriegel.  Causal transformer for estimating counterfactual outcomes.  _ArXiv_ , abs/2204.07258, 2022.
  * Miao et al. (2016) W. Miao, Z. Geng, and E. T. T. Tchetgen.  Identifying causal effects with proxy variables of an unmeasured confounder.  _Biometrika_ , 105 4:987–993, 2016.  doi: 10.1093/BIOMET/ASY038.  URL https://www.semanticscholar.org/paper/09341c95b87b2ab7c9a8a063b786f2c396687a10.
  * Miladinović et al. (2019) Đ. Miladinović, M. W. Gondal, B. Schölkopf, J. M. Buhmann, and S. Bauer.  Disentangled state space representations.  _arXiv preprint arXiv:1906.03255_ , 2019.
  * Morgan (2013) S. L. Morgan.  _Handbook of causal analysis for social research_.  Springer, 2013.
  * Mueller and Pearl (2023) S. Mueller and J. Pearl.  Personalized decision making–a conceptual introduction.  _Journal of Causal Inference_ , 11(1):20220050, 2023.
  * Müller (1997) A. Müller.  Integral probability metrics and their generating classes of functions.  _Advances in applied probability_ , 29(2):429–443, 1997.
  * Nie and Wager (2021) X. Nie and S. Wager.  Quasi-oracle estimation of heterogeneous treatment effects.  _Biometrika_ , 108(2):299–319, 2021.
  * Owen (2001) A. B. Owen.  _Empirical likelihood_.  Chapman and Hall/CRC, 2001.
  * Paszke et al. (2019) A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala.  Pytorch: An imperative style, high-performance deep learning library.  In _Neural Information Processing Systems_ , 2019.
  * Pearl (2009) J. Pearl.  _Causality_.  Cambridge university press, 2009.
  * Peters et al. (2017) J. Peters, D. Janzing, and B. Schölkopf.  _Elements of causal inference: foundations and learning algorithms_.  The MIT Press, 2017.
  * Potter et al. (2004) S. M. Potter, P. M. Robinson, S. Sakata, A. Smith, and D. L. Thornton.  Estimation and inference of impulse responses by local projections.  2004\.
  * Ranganath and Perotte (2018) R. Ranganath and A. Perotte.  Multiple causal inference with latent confounding.  _arXiv preprint arXiv:1805.08273_ , 2018.
  * Robins (1997) J. M. Robins.  Causal inference from complex longitudinal data.  In _Latent variable modeling and applications to causality_ , pages 69–117. Springer, 1997.
  * Robins and Hernán (2009a) J. M. Robins and M. A. Hernán.  Estimation of the causal effects of time-varying exposures.  _Longitudinal data analysis_ , 553:599, 2009a.
  * Robins and Hernán (2009b) J. M. Robins and M. A. Hernán.  Estimation of the causal effects of time-varying exposures.  _Longitudinal data analysis_ , 553:599, 2009b.
  * Robins et al. (2000) J. M. Robins, M. A. Hernan, and B. Brumback.  Marginal structural models and causal inference in epidemiology, 2000.
  * Rosenbaum and Rubin (1983) P. R. Rosenbaum and D. B. Rubin.  The central role of the propensity score in observational studies for causal effects.  _Biometrika_ , 70:41–55, 1983.
  * Rubin (2005) D. B. Rubin.  Causal inference using potential outcomes: Design, modeling, decisions.  _Journal of the American Statistical Association_ , 100(469):322–331, 2005.
  * Schuler et al. (2018) A. Schuler, M. Baiocchi, R. Tibshirani, and N. Shah.  A comparison of methods for model selection when estimating individual treatment effects.  _arXiv preprint arXiv:1804.05146_ , 2018.
  * Shalit (2020) U. Shalit.  Can we learn individual-level treatment policies from clinical data?  _Biostatistics_ , 21(2):359–362, 2020.
  * Shalit et al. (2017) U. Shalit, F. D. Johansson, and D. Sontag.  Estimating individual treatment effect: generalization bounds and algorithms.  In _International Conference on Machine Learning_ , pages 3076–3085. PMLR, 2017.
  * Shi et al. (2019) C. Shi, D. Blei, and V. Veitch.  Adapting neural networks for the estimation of treatment effects.  _Advances in neural information processing systems_ , 32, 2019.
  * Simchoni and Rosset (2021) G. Simchoni and S. Rosset.  Using random effects to account for high-cardinality categorical features and repeated measures in deep neural networks.  _Advances in Neural Information Processing Systems_ , 34, 2021.
  * Sinkhorn (1967) R. Sinkhorn.  Diagonal equivalence to matrices with prescribed row and column sums.  _The American Mathematical Monthly_ , 74(4):402–405, 1967.
  * Spirtes (2010) P. Spirtes.  Introduction to causal inference.  _Journal of Machine Learning Research_ , 11(5), 2010.
  * Sriperumbudur et al. (2009) B. K. Sriperumbudur, K. Fukumizu, A. Gretton, B. Schölkopf, and G. R. Lanckriet.  On integral probability metrics,$\backslash$phi-divergences and binary classification.  _arXiv preprint arXiv:0901.2698_ , 2009.
  * Sutskever et al. (2014) I. Sutskever, O. Vinyals, and Q. V. Le.  Sequence to sequence learning with neural networks.  _Advances in neural information processing systems_ , 27, 2014.
  * Tan (2010) Z. Tan.  Bounded, efficient and doubly robust estimation with inverse weighting.  _Biometrika_ , 97(3):661–682, 2010.
  * Tran et al. (2020) M.-N. Tran, N. Nguyen, D. Nott, and R. Kohn.  Bayesian deep net glm and glmm.  _Journal of Computational and Graphical Statistics_ , 29(1):97–113, 2020.
  * VanderWeele and Robins (2007) T. J. VanderWeele and J. M. Robins.  Four types of effect modification: A classification based on directed acyclic graphs.  _Epidemiology_ , 18:561–568, 2007.
  * Vansteelandt and Joffe (2014) S. Vansteelandt and M. Joffe.  Structural nested models and g-estimation: the partially realized promise.  2014\.
  * Villani et al. (2009) C. Villani et al.  _Optimal transport: old and new_ , volume 338.  Springer, 2009.
  * Viviano and Bradic (2021) D. Viviano and J. Bradic.  Dynamic covariate balancing: estimating treatment effects over time.  2021\.
  * Wager and Athey (2018) S. Wager and S. Athey.  Estimation and inference of heterogeneous treatment effects using random forests.  _Journal of the American Statistical Association_ , 113:1228 – 1242, 2018.
  * Wang et al. (2022) X. Wang, S. Lyu, X. Wu, T. Wu, and H. Chen.  Generalization bounds for estimating causal effects of continuous treatments.  _Advances in Neural Information Processing Systems_ , 35:8605–8617, 2022.
  * Wang and Blei (2019) Y. Wang and D. M. Blei.  The blessings of multiple causes.  _Journal of the American Statistical Association_ , 114(528):1574–1596, 2019.
  * Wilcoxon (1992) F. Wilcoxon.  Individual comparisons by ranking methods.  In _Breakthroughs in Statistics: Methodology and Distribution_ , pages 196–202. Springer, 1992.
  * Xiong et al. (2019) Y. Xiong, H. J. Kim, and V. Singh.  Mixed effects neural networks (menets) with applications to gaze estimation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 7743–7752, 2019.
  * Zhang et al. (2020a) Y. Zhang, A. Bellot, and M. Schaar.  Learning overlapping representations for the estimation of individualized treatment effects.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1005–1014. PMLR, 2020a.
  * Zhang et al. (2020b) Y. Zhang, A. Bellot, and M. Schaar.  Learning overlapping representations for the estimation of individualized treatment effects.  In _International Conference on Artificial Intelligence and Statistics_ , pages 1005–1014. PMLR, 2020b.
  * Zhou and Wodtke (2018) X. Zhou and G. T. Wodtke.  Residual balancing weights for marginal structural models: with application to analyses of time-varying treatments and causal mediation.  _arXiv: Applications_ , 2018.
  * Zubizarreta (2015) J. R. Zubizarreta.  Stable weights that balance covariates for estimation with incomplete outcome data.  _Journal of the American Statistical Association_ , 110(511):910–922, 2015.


# Performance of Polarization-based Stereoscopy Screens

Xiaozhu Zhang  Kristian Hantke  Max Planck Institute for Dynamics and Self-
Organization (MPIDS), 37077 Goettingen, Germany  Cornelius Fischer  Georg-
August-Universität Göttingen, 37077 Göttingen, Germany  Matthias Schröter
<EMAIL_ADDRESS>Max Planck Institute for Dynamics and Self-
Organization (MPIDS), 37077 Goettingen, Germany

###### Abstract

The screen is a key part of stereoscopic display systems using polarization to
separate the different channels for each eye. The system crosstalk,
characterizing the imperfection of the screen in terms of preserving the
polarization of the incoming signal, and the scattering rate, characterizing
the ability of the screen to deliver the incoming light to the viewers,
determine the image quality of the system. Both values will depend on the
viewing angle. In this work we measure the performance of three silver screens
and three rear-projection screens. Additionally, we measure the surface
texture of the screens using white-light interferometry. While part of our
optical results can be explained by the surface roughness, more work is needed
to understand the optical properties of the screens from a microscopic model.

stereoscopic projection and polarization and screen and system crosstalk and
ghosting and surface texture

##  I Introduction

Displaying 3D content is not only an important issue in the entertainment
industry, it is also of increasing importance in science where new numeric and
experimental methods have created a wealth of three-dimensional datasets. Many
stereoscopic display system are based on polarization filtering: the visual
information for each eye is oppositely polarized, projected to and scattered
or transmitted by the screen, and finally filtered by the viewer’s glasses
which consist of two polarizers admitting only the correctly polarized light
to each eye iizuka:06 ; janssen:08 ; kim:10 . There are two options for
polarization filtering: linear and circular polarized light. While linear
polarizers are simpler to manufacture, circular polarization has the advantage
that head tilting will not impair the quality of the image

An ideal screen would completely preserve the polarization of the incoming
light. However, in practice there is always some amount of “ghosting”
resulting from the change of polarization at the screen. A measure for
ghosting is the system crosstalk $C$. It is defined as the ratio between the
intensity of light that leaks from the unintended channel to the intended one
and the intensity of the intended channel woods:11 . According to measurements
of Huang et al. huang:03 , the maximal acceptable system crosstalk for a
typical viewer to still experience a stereo sensation is 0.1. (Lower values
down to $10^{-4}$ can still be detected by careful visual inspection). While
it is also known on a theoretical basis that the viewing angle will influence
the amount of system crosstalk richards:10 , to our knowledge no measurements
of the angle-dependent system crosstalk of different screen types have been
published up to now. Neither has the question been studied how the inclination
angle (between the incoming light from the projector and the surface normal of
the screen) influences the system crosstalk.

A second measure for the quality of a screen is the brightness of the image,
which depends on the amount and angular distribution of the reflectance (for
silver screens) or transparency (for rear-projection screens) of the screen.
For silver screens this is typically quoted as the screen gain, the intensity
measured at normal incidence normalized by the intensity of a Lambertian
source brennesholtz:08 . Here we measure the angle dependent scattering rate
$S$ for both silver screens and rear-projection screens. $S$ is defined as the
ratio of the intensity received by a viewer in a certain angle to the
intensity of the incoming light, normalized by the solid angle.

In this paper we present measurements of the angular dependence of system
crosstalk and scattering rate for three samples of silver screens (labeled SS1
to SS3) and three rear-projection screens (RP1 to RP3). Additionally, we
determine the surface texture of the samples using white-light interferometry;
this information provides some qualitative insight into our optical results.

##  II Experimental Setup

|  |
---|---|---
(a) | (b) | (c)
Figure 1: Sketch of the experimental setup using circular polarized light:
panels a) and b) are used for silver screens and a) and c) for rear-projection
screens. The generation and the detection of linearly polarized light is
achieved by removing the devices colored in red.

Figure 1 shows the experimental setup used for measuring the angular
dependence of the system crosstalk and the scattering rate. Diode pumped solid
state lasers (DPGL-2050 from Photop and Verdi V5 SF from Coherent) with a
wavelength of 532nm were used as light sources for the experiments. Passing a
beam expander, the diameter of the laser beam was increased to 3.4mm (FWHM),
whereas the typical size of structural inhomogeneities on the screen surface
is at most a few hundred micrometers as shown below. This ensured that the
measured data for different spots on the screen are reproducible within $\pm
5$%.

The laser light was linearly polarized by passing a polarizer or circularly
polarized by passing an additional Babinet-Soleil compensator (from B. Halle).
The screen sample is irradiated by the laser at normal incidence and the
scattered laser light of the silver and rear projection screens is detected by
an detection unit in reflection (1b) and in transmission (1c), respectively.
The detection unit consists of a power meter (PM100D with sensor S130C from
Thorlabs), an analyzer and in case of the circular polarization an additional
quarter-wave plate (both from B. Halle). A long and narrow tube was placed in
front of the power meter, ensuring that only the photons are detected that
scatter from the irradiated spot on the screen along the viewing axis of the
sensor in a solid angle $\Omega$ of 2$\cdot$10-4 sr. The detection unit was
placed on a rotatable rail with the rotational axis being fixed in such a way
that the normal viewing axis of the sensor intercepts always with the
illuminated area on the screen during rotation. The viewing angle $\theta$ can
be varied from -20° to 80°. For silver screens the range of $\pm$6° is
inaccessible in order to not block the incoming beam.

In both cases of linear and circular polarization the incoming laser intensity
$I_{in}$ was measured just in front of the sample. Furthermore, the intensity
of the scattered light $I_{out}$ was measured for the intended channel with
the polarization being the same direction as the incoming one ($I_{out}^{sp}$,
analyzer and polarizer parallel) and for the unintended channel with the
polarization being the opposite direction ($I_{out}^{op}$, analyzer and
polarizer perpendicular) for different viewing angles $\theta$. From this data
one can compute the crosstalk $C(\theta)$:

$C(\theta)=\frac{I_{out}^{op}(\theta)}{I_{out}^{sp}(\theta)}$ (1)

and the scattering rate $S(\theta)$:

$S(\theta)=\frac{I_{out}^{sp}(\theta)+I_{out}^{op}(\theta)}{I_{in}\;\Omega}$
(2)

The precision of the measurement for the crosstalk depends strongly on the
purity of the initial laser polarization, whereas the scattering rate is not
affected within our measurement precision. Analyzing the crosstalk without any
screen sample (i.e. putting the laser directly in front of the analyzer
system) we found the lower resolution limit in the linear case to be less than
$3\times 10^{-3}$. In the circular case the degree of polarization results in
a lower resolution limit of $7.5\times 10^{-3}$.

The screen samples were obtained from the company Screenlab (Elmshorn,
Germany), their specifications and brand names are listed in table 1.

Table 1: Sample labels, brand names, manufacturer information on gain and transmission, and surface properties measured by white light interferometry: the root mean square roughness $R_{q}$ and the ratio between the surface area and the projected area $F$. Sample | brand name | gain | transmission | $R_{q}[\mu m]$ | $F$
---|---|---|---|---|---
RP1 | BS XRP3 |  | 41.8 | 3.8 | 1.03
RP2 | WS XRP3 |  | 88.8 | 3 | 1.13
RP3 | BS RP2 |  | 41.2 | 4.2 | 1.2
SS1 | SH120 | 2.4 |  | 5 | 2.3
SS2 | SF120 | 2.4 |  | 8 | 2.6
SS3 | WA160 | 1.3 |  | 22 | 2.2
Figure 2: Angular dependence of the crosstalk caused by the screen using a)
circular polarized light and b) linear polarization. The grey dash line
corresponds to the threshold for still acceptable stereo fusion according to
Huang et al. huang:03 .

The surface topography of the screens was measured using a ZeMapper whitelight
vertical scanning interferometer (Zemetrics, Tucson, USA): the focal plane of
an interference pattern is vertically scanned through the sample topography,
then a height map is calculated from the collected amplitude maps of the
interference patterns. The vertical resolution of the instrument is better
than 1 nm; the maximum field of view applied in this study is 1.4 mm. For more
information on the instrument see darbha:10 . Prior to the measurement the
rear projection screens where sputter coated with a 40 nm gold layer to
increase surface reflectivity.

##  III Results

Figure 2 displays the system crosstalk of the six screen samples, both with
circular and linear polarized light. Based on the criterion found by Huang et
al. huang:03 , all screens allow stereo vision for viewing angles $\theta$
smaller than 40° (circular) or 48° (linear). In practice this range will be
smaller due to the additional crosstalk originating from the glasses and the
inclination angle of the incoming light; the latter effect will be described
below.

In general, the screens seem to fall into two categories; they are either
optimized for a large range of acceptable crosstalk or a minimized crosstalk
at small $\theta$. In both categories the silver screens are outperformed by
the rear projection screens: RP3 has a smaller $C$ at small $\theta$ than SS1
while RP1 has a broader range of acceptable viewing angles than SS3.

Regarding the polarization mode, linear polarization has for each screen a
clear advantage over circular. $C_{circ}/C_{lin}$ measured at $\theta$ = 10°
varies between 1.1 (RP1) and 4 (RP3) as shown in figure 4. Please observe,
that our measurements of the crosstalk of RP3 at small angles might be limited
by our experimental resolution.

Figure 3: Ratio of crosstalk measured with circular and linear polarization of
the incoming light.

Figure 4: Dependence of the crosstalk on the inclination angle $\phi$.
Measured on screen SS1 using linear polarized light.

Under real world conditions it is quite likely that the incoming light itself
will have an inclination angle $\phi$ to the surface normal of the screen. To
quantify the additional crosstalk created this way, we modified the
experimental setup by adding a periscope in front of the polarizer. Figure 4
shows the crosstalk for sample SS1 in the case of linear polarization for
inclination angles between 0°and 15°. While there is a clear increase of
crosstalk with $\phi$, the range of acceptable viewing angles $\theta$ is
reduced by only 5°.

Figure 5: Angular dependence of the scattering rate measured with circular
polarization. Scattering rate values for linear polarization agree within 5.4
percent.

The scattering rates $S$ of the screen samples with circular polarization are
shown in figure 5. Deviations of $S$ measured with linear or the circular
polarization are within our errorbars. For high luminosities at small viewing
angles SS1 and RP2 are the best choice, in terms of best homogeneity SS3 comes
closest to a Lambertian source.

Figure 6: Perspective images of the surface texture of the screen samples.
Please note the different horizontal and vertical scales for sample SS3.
Images contain between 0.7% (RP1) and 30% (SS3) interpolated pixels.

From a theoretical side the performance of a screen will depend both on its
material and its surface texture jin:10 . While we do not have information on
the electromagnetic properties of the screen material, the surface texture can
be measured with white light interferometry. Perspective images of the
surfaces of SS1, SS3, RP1, RP3 are shown in figure 6. The RMS (root mean
square) roughness $R_{q}$ and the ratio between the surface area and the
projected area $F$ of all six screen samples is listed in table 1.

A comparison of the angular dependence of $S$ with these values hints at
$R_{q}$ as a predictor for the deviation from a Lambertian source. This is
particularly shown by SS3 which has by far the highest value of $R_{q}$ and
the smallest $\theta$ dependence of $S$. Regarding the system crosstalk a
similar correlation between $R_{q}$ and the slope of $C$ at large angles might
exist. On the other side we do not find a clear correlation between the
optical properties and $F$.

##  IV Conclusion

All screens allow effective stereo projection for viewing angles up to 40°. At
larger angles the crosstalk of rear projection screens is considerably smaller
than that of silver screens. Also for each screen the crosstalk was larger
with circular polarization than with linear. However, when planing a display
system additional factors have to be taken into account like the available
space behind the screen or the sensitivity of the system against the viewers
tilting their heads. Consequentially, no optimal solution for all possible
scenarios exist. While the roughness of the screens influences their large
viewing angle behavior, clearly more research is needed for a quantitative
understanding.

Acknowledgements: We would like to thank Günter Daszinnies from the company
Screenlabs for providing the test samples.

## References

  * (1) Brennesholtz, M.S., Stupp, E.H.: Projection Displays.  Wiley (2008)
  * (2) Darbha, G., Schäfer, T., Heberling, F., Lüttge, A., Fischer, C.: Retention of latex colloids on calcite as a function of surface roughness and topography.  Langmuir 26, 4743–4752 (2010)
  * (3) Huang, K.C., Yuan, J.C., Tsai, C.H., Hsueh, W.J., Wang, N.Y.: A study of how crosstalk affects stereopsis in stereoscopic displays.  In: Proceedings of SPIE-IS&T, vol. 5006, pp. 247–253 (2003)
  * (4) Iizuka, K.: Welcome to the wonderful world of 3D: introduction, principles and history.  Optics and Photonics News 17, 42–51 (2006)
  * (5) Janssen, J.K.: 3D 2.0, Neuer Anlauf für Stereoskopie im Kino.  c’t 16, 72–75 (2008)
  * (6) Jin, L., Kasahara, M., Gelloz, B., Takizawa, K.: Polarization properties of scattered light from macrorough surfaces.  Optics Letters 35, 595–597 (2010)
  * (7) Kim, S.C., Kim, E.S.: Performance analysis of stereoscopic three-dimensional projection display systems.  3D Research 1, 1–16 (2010)
  * (8) Richards, M., Schnuelle, D.: The effective gain of a projection screen in an auditorium.  SMPTE Motion Imaging Journal 119, 62 –67 (2010)
  * (9) Woods, A.J.: How are crosstalk and ghosting defined in the stereoscopic literature?  In: Proceedings of SPIE-IS&T, vol. 7863, p. 78630Z (2011)


††thanks: These authors contributed equally to this work.††thanks: These
authors contributed equally to this work.

# Tuning Colloidal Reactions

Ryan Krueger  School of Engineering and Applied Sciences, Harvard University,
29 Oxford St, Cambridge, MA 02138, USA  Ella King  Department of Physics, New
York University, 726 Broadway, New York, NY 10003, USA  Simons Junior Fellow,
160 5th Ave, New York, NY 10010, USA  Michael Brenner
<EMAIL_ADDRESS>School of Engineering and Applied Sciences, Harvard
University, 29 Oxford St, Cambridge, MA 02138, USA

###### Abstract

The precise control of complex reactions is critical for biological processes
ranging from cell division to metabolism. Synthetic analogues of living
materials suffer from our inability to tune chemical reactions with precise
outcomes. Here, we leverage differentiable simulators to design nontrivial
reaction pathways in colloidal assemblies. By optimizing interactions between
reactants and substrates, we achieve controlled disassembly of octahedral and
icosahedral shells. As a potential engineering target, we design a reaction
that provokes the release of a small particle trapped in a shell.

††preprint: APS/123-QED

Both living and non-living physical systems exhibit complex, dynamical
behavior, ranging from repair, locomotion, and catalysis. Fundamentally, these
complex dynamical behaviors arise from sequences of reactions, in which a set
of substances (i.e. the reactants) are transformed into a set of different
substances (i.e. the products).

A rich body of theoretical work aims to characterize and tune systems of
interacting agents spanning a range of system descriptions [1, 2, 3]. There
have been exciting advances in creating tunable reactions that are
experimentally realizable [4, 5, 6, 7, 8, 9]. However, the models that provide
an experimentally-relevant level of detail are typically bespoke and
application-specific, necessitating entirely new research programs for each
new reaction. Conversely, general theoretical models have provided deep
insights into the basis of self-limiting growth [10, 11], self-replication
[12],catalysis [13, 14], and much more [15, 16, 17, 18, 19]. These models,
however, are largely too abstract to inform experimental design. Moreover,
many reactions that yield desirable dynamical behavior are incredibly finely-
tuned, making engineering these reactions untenable with most traditional
design approaches. To design precise reactions, it is necessary to carry out
inverse design, whereby one designs components and their interactions to
achieve a given reaction. While inverse design has successfully applied to
self assembly [20, 21, 22, 23, 24], its application to reaction pathways is
much more difficult owing to the need to choose design parameters to favor
particular dynamical trajectories. The advent of differentiable simulators
[24], powered by software libraries developed for machine learning [25], has
opened up the possibilty of directly designing reactions.

Here, we use _differentiable molecular dynamics (MD)_ to carry out the inverse
design of complex reactions. The differentiability makes it possible to
compute analytical gradients with respect to trajectories during the
interaction of components and perform gradient-based optimization with respect
to parameters to achieve a target reaction. To demonstrate the ability
framework to design nontrivial reactions, we consider the controlled
disassembly of colloidal structures, whereby one particle is extracated from
an otherwise complete shell of colloidal particles. Disassembly is central to
the dynamic functions of living systems, such as defect repair, self-
replication, and catalysis. Most existing examples of controlled disassembly
in synthetic systems rely on external forcing to drive the disassembly process
[26, 27, 28, 29], which provides a direct pathway to tuning behavior. However,
for many engineering applications including those inherent to living systems,
the use of external fields is limiting. Alternatively, living systems
typically rely on local energy consumption (e.g. biological enzymes consuming
ATP) rather than global fields, but the synthetic design of these systems is
significantly more complex.

In this paper, we design for the controlled disassembly of both an octahedral
shell and an icosahedral shell. Importantly, our disassembly mechanism does
not rely on external forcing and occurs without consuming energy. As a model
for potential engineering applications, we apply our disassembly mechanism to
release a small particle trapped inside an icosahedral shell. Controlled
disassembly serves as a striking example of a complex reaction because the
reaction requires a finely-tuned interaction energy to keep the remaining
shell stabilized while performing the desired disassembly. Moreover, we
consider a fixed substrate and only parameterize an external structure that
acts on this fixed substrate, enabling control over disassembly without
modifying critical components of the target reaction.

Figure 1:  An overview of our method for tuning the interaction potential of
an spider to achieve a desired reaction. A. The desired reaction, in which a
target vertex is removed from an assembled octahedron. B.  A candidate
mechanism for this reaction, in which an external structure called the
“spider” extracts the target vertex via an attractive potential and detaches
from the remaining shell. C.  The parameterization of the geometry of the
spider and its interaction potential with the shell. The red particle is the
“head” particle, situated above the four black “base” particles that
constitute the ring. The interaction energy between the spider and the shell
is depicted as a green triangle. We optimize over all labeled parameters, as
well as the cutoff of the interaction energy which is not depicted. D.  A
high-level depiction of our optimization pipeline, in which analytic gradients
are computed via a differentiable molecular dynamics simulator and parameters
of the spider are updated via gradient descent.  Figure 2:  An overview of the
optimizations for the octahedral shell, performed in the limit of weak and
strong initial interactions between the shell and spider (panels A and B,
respectively). Insets depict representative states after 10,000 MD steps with
the parameters at the corresponding optimization iteration. A. In the limit of
weak initial interaction between the shell and spider, the initial spider
simply diffuses away from the shell. By the 500th iteration, the spider
geometry and interaction energy are optimized to extract the target vertex
while still diffusing away from the remaining vertices and leaving them
intact. B.  In the limit of strong initial interaction between the shell and
spider, the initial spider extracts the target vertex but does not diffuse
away and disturbs the integrity of the remaining shell. As optimization
progresses, the interaction is tuned to only extract the target vertex while
not disrupting the remaining vertices. Upon convergence, the spider geometry
and interaction energy are tuned to maintain extraction while diffusing away
from the remaining shell. When plotting the individual loss terms, we plot the
logarithm of the remaining energy term rather than its raw value such that the
scales of the two terms are comparable for the entire range of the figure.

Results.—We implement controlled disassembly example in a patchy particle
system. Patchy particles have long been used to emulate interactions in soft
materials[30, 31], and offer tremendous flexibility of designed interactions,
especially when patches can be programmed to interact with each other with
designed strength and range. Optimizing patchy particle systems to achieve
specific behaviors has been made possible by the recent development of patchy
particle simulations within a differentiable library [24].The ability to carry
out this type of optimization within a molecular dynamics simulations opens up
enormous possibilities for the design of novel reactions. In particular, we
aim to remove a single particle from a shell composed of patchy particles
without disrupting the remaining shell structure. We tune disassembly without
changing any properties of the shell itself, by introducing an external
reactant that interacts with the shell to disassemble it in the desired
manner. We term the external reactant a “spider” due to the structure of the
object.

The shell is modelled as a collection of patchy particles where each patch
corresponds to a contact with a neighboring particle. Patches interact via a
Morse potential ($\epsilon_{V}=10.0$, $\alpha_{V}=5.0$) and vertex centers
interact via soft-sphere repulsion ($\epsilon_{ss}=10,000$). Importantly, the
geometry and interaction energy of a shell vertex are fixed throughout an
optimization: only features of the spider are optimized.

The spider is modeled as a rigid body composed of a ring of “base” particles
and a “head” particle that sits above the ring along its symmetry axis. The
base particles are connected by repulsive bars and the head is connected to
the base particles particles by the same repulsive bars, making the entire
structure a cage-like object that is open on one end. The head particle
interacts with shell vertices via a Morse potential whereas base particles and
connecting bars interact with shell vertices via soft sphere repulsion. Unlike
the shell, the geometry and interaction energy of the spider are parameters of
the optimization. See Figure 1C for an overview of this parameterization. All
the interaction energies in our system are parameterized with simple, physics-
based potentials.

We optimize over 8 parameters that characterize the geometry of the spider and
its interaction with the shell via differentiable MD (see Figure 1D). Our loss
term is constructed via two competing terms: one that rewards a final state in
which the target vertex is extracted, and one that penalizes a strong
interaction between the spider and non-target vertices. Consider a shell
comprised of a collection of vertices
$V=\\{\overrightarrow{v}_{1},\overrightarrow{v}_{2},\cdots,\overrightarrow{v}_{n}\\}$
where $n=8$ and $n=12$ for the octahedron and icosahedron, respectively. We
seek to extract a target vertex $\overrightarrow{v}_{j}$ from the shell while
leaving the remaining shell $V\setminus\overrightarrow{v}_{j}$ intact. We can
measure the degree to which $\overrightarrow{v}_{j}$ is successfully extracted
via the following expression:

$\displaystyle\mathcal{L}_{\text{extract}}(V)=-\sum_{i\neq
j}d(\overrightarrow{v}_{i},\overrightarrow{v}_{j})$ (1)

where $d(\overrightarrow{v}_{i},\overrightarrow{v}_{j})$ denotes the Euclidean
distance between vertices $\overrightarrow{v}_{i}$ and
$\overrightarrow{v}_{j}$. Note that the negative sign as we formulate our
optimization problem to minimize the loss. Secondly, we can minimize the
interaction energy between the spider and the non-target vertices of the
shell,

$\displaystyle\mathcal{L}_{\text{remain}}(V,\overrightarrow{h})=\left(\sum_{i\neq
j}U_{\text{ss}}(\overrightarrow{h},\overrightarrow{v}_{i})\right)^{2}$ (2)

where $U_{\text{ss}}(\overrightarrow{h},\overrightarrow{v}_{i})$ represents
the interaction energy between the spider head, $\overrightarrow{h}$, and a
shell vertex $\overrightarrow{v}_{i}\in V$ as depicted in Figure 1C. For the
total loss, we evaluate the interaction between the spider and the non-target
shell vertices with respect to the initial state while we calculate the
spider’s extraction of the target vertex with respect to the final state. We
evaluate the remaining energy term on the initial state because this term is
only well defined when evaluated on the same configuration at each iteration.
In all simulations, the spider is initialized bound to the target vertex and
the system is integrated for 1000 timesteps (see Appendix).

Figure 3:  An overview of our optimization procedure to achieve controlled
release of a particle from an icosahedral shell. A. The target reaction, in
which a spider extracts a target vertex and the cargo is released from the
remaining stable shell. B. The loss values for two optimizations initialized
with weak ($log(\epsilon_{H})=5.5$, $\alpha_{H}=1.5$) and strong
($log(\epsilon_{H})=9.2$, $\alpha_{H}=1.5$) interactions between the spider
and shell. C. Representative trajectories of the spider-vertex complex and
cargo for a simulation with 10,000 timesteps. The cargo remains within its
container until the spider-vertex complex diffuses from the remaining shell,
upon which the cargo can exit through the top.

Octahedron.—We first demonstrate our method on an octahedron, a simple but
nontrivial platonic solid. We explore two limits of our optimization
procedure: weak and strong initial shell-spider interactions. Figure 2
summarizes our results.

First, we perform an optimization where the spider is initialized to interact
weakly with the shell vertices ($log(\epsilon_{H})=3.5$, $\alpha_{H}=1.5$).
Under this weak initial interaction, the spider simply diffuses away from the
shell at long timescales without extracting the target vertex. Then,
throughout the course of the optimization, we observe variable changes
consistent with increasing the interaction between the spider head and the
shell – $\epsilon_{H}$ increases, the head height decreases, and the head
particle radius increases. However, we also observe less intuitive changes –
$\alpha_{H}$ increases (decreasing the range of the Morse potential), the
shell base radius decreases, and the base particle radius increases. This
suggests tightly coupled, nontrivial parameter changes drive the extraction of
the target vertex while maintaining minimal interaction with the remaining
shell.

Next, we performed an optimization in the opposite limit in which the spider
is initialized to interact strongly with the shell ($log(\epsilon_{H})=8.5$,
$\alpha_{H}=1.0$). Initially, this interaction is so strong that the spider
not only extracts the target vertex but it also disrupts the remaining shell.
This disruption can be seen in the large value of the remaining energy loss
term, which penalizes the energy between the spider head and non-target
vertices. Throughout the optimization, we observe variable changes consistent
with tuning the interaction strength to maintain extraction while minimizing
off-target interactions – $\epsilon_{H}$ decreases, $\alpha_{H}$ increases,
the head radius decreases, the head height increases, and the base particle
radius increases. When evaluated on longer simulations, the converged
parameter set also achieves diffusion of the spider-vertex complex from the
remaining, undisturbed shell.

While only single optimizations are presented in Figures 2 and 3, we
additionally performed ensembles of optimizations over different initial seeds
for the random number generator and over different random perturbations to the
initial $log(\epsilon_{H})$ parameter in both the high and low energy regimes.
Despite the wide variance in initial parameters, we observe low variance in
the converged parameter values in all cases (see Figure 4). This may suggest
that the solution set is narrow.

Contrasting these two optimization regimes reveals the inherent delicacy in
tuning the spider to achieve extraction of a target vertex and subsequent
diffusion of the vertex away from the shell. The spider-shell interaction must
be sufficiently strong to extract the target vertex, but be simultaneously
weak enough to not disturb the non-target vertices and to diffuse away from
the shell within the timescale of our simulations. This tension can be seen in
the behavior of the loss terms in each optimization. In the weak-interaction
limit, the term penalizing interactions with non-target vertices remains
negligible while the extraction term drives optimization; in the strong-
interaction limit it is the same energy-penalizing term that dominates the
loss.

Figure 4:  An evaluation of the variation in the optimized parameters for each
optimization regime for the octahedral shell. For both the low- and high-
energy initial parameter regimes, two sets of five optimizations were
performed either with different random keys or with perturbations to the
initial value of $\log(\epsilon_{H})$ drawn from $\mathcal{N}(0,0.1)$.
Optimization was performed for 500 and 1250 iterations for the low- and high-
energy limits, respectively. Figures A-D depict box plots of the final
parameters across each set of 5 optimizations. For both different random keys
and perturbed initial $\log(\epsilon_{H})$, we observe low variance in the
converged parameter values. The lower variance in the high-energy limit in
comparison to the low-energy limit may result from the higher number of
iterations.

Icosahedron and Controlled Release.—We now increase in complexity to a model
of a biologically relevant function: controlled release of a particle from an
icosahedral shell. We initialize our system with a small particle trapped
inside a fully formed icosahedral shell. We then optimize a similar spider
structure to perform controlled disassembly of the icosahedral shell,
resulting in the controlled release of the small particle. This serves as a
toy example of a potential target for engineering applications, such as drug
delivery via a viral shell.

To match the symmetry of the icosahedron, we use a spider with five legs
rather than four. Otherwise, optimization proceeds in the same fashion as for
the octahedral shell. We consider the same two limits: one with high
interaction strength between the spider and the shell and one with low
interaction strength between the spider and the shell (see Figure 3). The same
delicate tradeoff between extraction and diffusion manifests with the
icosahedron, and our automatic differentiation-based optimization deftly
navigates this tradeoff.

Discussion.—In this Letter, we have demonstrated the design of nontrivial
reactions via differentiable simulation. We considered the case of controlled
disassembly, in which there is an inherent tension between initiating
disassembly and maintaining the integrity of the remaining substructure. We
first consider the controlled disassembly of an octahedral shell composed of
patchy particles and show how the parameters governing an external structure
(i.e. the “spider”) can be finely tuned to minimize a loss function
representing this inherent tension. We then increase in complexity,
demonstrating the controlled release of a small particle from an icosahedral
shell.

By optimizing directly with respect to the numerically integrated dynamics,
our method is general enough to study a wide range of systems. Foremost, our
method may enable the experimental realization of theoretical models that were
otherwise limited by an inability to finely tune interaction energies. For
example, Ref. [12] introduces a model of self-replicating colloidal clusters
in which kinetic traps can be avoided by tuning the interaction energies, but
dissociation of a new cluster from its parent (a necessary step for
replication) required an artificial trigger event in numerical simulations. In
contrast, our designed parameters lead to spontaneous dissociation of the
spider-vertex complex away from the remaining shell. The computational
flexibility of the method could also easily allow restricting the parameter
regime to experimentally realizable interactions. This could be done for DNA
coated colloids, for instance, by optimizing the DNA sequences that define the
interaction strength [32, 33].

A number of numerical instabilities can arise while optimizing for the types
of reactions we consider. The primary limitation we observe is that gradients
become unstable and very large for long simulations. There are several
possible approaches to reducing instability in gradients of long simulations.
One standard method to mitigate such instabilities in the context of
differentiable programming is gradient clipping [34, 35]. One could also
decrease the total number of timesteps by training an emulator to resolve the
dynamics with a larger timestep than is otherwise possible with standard
integrators, following similar work for deterministic systems [36, 37, 38]. An
alternative approach to decreasing the number of timesteps would be to
integrate differentiable simulations with enhanced sampling methods. As a
result, one could sample low probability events without the need for long
simulation times.

We rely on gradient-based optimization due to its scalability and performance.
The method scales naturally to larger systems since the calculation of a
gradient via automatic differentiation only requires a single simulation,
while at the same time reverse-mode automatic differentiation scales
efficiently with the number of parameters [39]. Calculating the analytical
gradient with respect to a simulation explicitly captures interdependency
between parameters, essential to efficiently tuning complex behavior.

The experimental realization of theoretical models of complex dynamical
behavior has long been limited by the ability to finely tune the
parameterization of such models. This method serves as a potential bridge
between these two worlds, permitting the delicate tradeoffs inherent in
parameterizations of complex reactions to be navigated efficiently using
physics-based potentials and dynamics.

###### Acknowledgements.

We thank Sam Schoenholz for his work developing JAX-MD and the members of the
Brenner Group at Harvard University for helpful discussions. This material is
based upon work supported by the Office of Naval Research (ONR
N00014-17-1-3029, ONR N00014-23-1-2654), the NSF Grant DMR-1921619, the NSF AI
Institute of Dynamic Systems (#2112085), and a grant from the Simons
Foundation (#1141499).

## References

  * van Anders  _et al._ [2015] G. van Anders, D. Klotsa, A. S. Karas, P. M. Dodd, and S. C. Glotzer, Digital alchemy for materials design: Colloids and beyond, Acs Nano 9, 9542 (2015).
  * Zeravcic  _et al._ [2017] Z. Zeravcic, V. N. Manoharan, and M. P. Brenner, Colloquium: Toward living matter with colloidal particles, Reviews of Modern Physics 89, 031001 (2017).
  * Sherman  _et al._ [2020] Z. M. Sherman, M. P. Howard, B. A. Lindquist, R. B. Jadrich, and T. M. Truskett, Inverse methods for design of soft materials, The Journal of chemical physics 152 (2020).
  * Truong  _et al._ [2022] V. X. Truong, L. L. Rodrigues, and C. Barner-Kowollik, Light-and mechanic field controlled dynamic soft matter materials, Polymer Chemistry 13, 4915 (2022).
  * Li  _et al._ [2019] Q. Li, A. P. Schenning, and T. J. Bunning, Light-responsive smart soft matter technologies (2019).
  * Wang  _et al._ [2012] Y. Wang, Y. Wang, D. R. Breed, V. N. Manoharan, L. Feng, A. D. Hollingsworth, M. Weck, and D. J. Pine, Colloids with valence and specific directional bonding, Nature 491, 51 (2012).
  * McMullen  _et al._ [2022] A. McMullen, M. Muñoz Basagoiti, Z. Zeravcic, and J. Brujic, Self-assembly of emulsion droplets through programmable folding, Nature 610, 502 (2022).
  * Niu  _et al._ [2019] R. Niu, C. X. Du, E. Esposito, J. Ng, M. P. Brenner, P. L. McEuen, and I. Cohen, Magnetic handshake materials as a scale-invariant platform for programmed self-assembly, Proceedings of the National Academy of Sciences 116, 24402 (2019).
  * Schade  _et al._ [2013] N. B. Schade, M. C. Holmes-Cerfon, E. R. Chen, D. Aronzon, J. W. Collins, J. A. Fan, F. Capasso, and V. N. Manoharan, Tetrahedral colloidal clusters from random parking of bidisperse spheres, Physical review letters 110, 148303 (2013).
  * Murugan  _et al._ [2015a] A. Murugan, J. Zou, and M. P. Brenner, Undesired usage and the robust self-assembly of heterogeneous structures, Nature communications 6, 6203 (2015a).
  * Zeravcic  _et al._ [2014] Z. Zeravcic, V. N. Manoharan, and M. P. Brenner, Size limits of self-assembled colloidal structures made using specific interactions, Proceedings of the National Academy of Sciences 111, 15918 (2014).
  * Zeravcic and Brenner [2014] Z. Zeravcic and M. P. Brenner, Self-replicating colloidal clusters, Proceedings of the National Academy of Sciences 111, 1748 (2014).
  * Zeravcic and Brenner [2017] Z. Zeravcic and M. P. Brenner, Spontaneous emergence of catalytic cycles with colloidal spheres, Proceedings of the National Academy of Sciences 114, 4342 (2017).
  * Muñoz-Basagoiti  _et al._ [2023] M. Muñoz-Basagoiti, O. Rivoire, and Z. Zeravcic, Computational design of a minimal catalyst using colloidal particles with programmable interactions, Soft Matter 19, 3933 (2023).
  * Keim  _et al._ [2019] N. C. Keim, J. D. Paulsen, Z. Zeravcic, S. Sastry, and S. R. Nagel, Memory formation in matter, Reviews of Modern Physics 91, 035002 (2019).
  * Murugan  _et al._ [2015b] A. Murugan, Z. Zeravcic, M. P. Brenner, and S. Leibler, Multifarious assembly mixtures: Systems allowing retrieval of diverse stored structures, Proceedings of the National Academy of Sciences 112, 54 (2015b).
  * Holmes-Cerfon  _et al._ [2013] M. Holmes-Cerfon, S. J. Gortler, and M. P. Brenner, A geometrical approach to computing free-energy landscapes from short-ranged potentials, Proceedings of the National Academy of Sciences 110, E5 (2013).
  * Holmes-Cerfon [2016] M. C. Holmes-Cerfon, Enumerating rigid sphere packings, siam REVIEW 58, 229 (2016).
  * Murugan and Vaikuntanathan [2017] A. Murugan and S. Vaikuntanathan, Topologically protected modes in non-equilibrium stochastic systems, Nature communications 8, 13881 (2017).
  * Torquato [2009] S. Torquato, Inverse optimization techniques for targeted self-assembly, Soft Matter 5, 1157 (2009).
  * Chen  _et al._ [2018] D. Chen, G. Zhang, and S. Torquato, Inverse design of colloidal crystals via optimized patchy interactions, The Journal of Physical Chemistry B 122, 8462 (2018).
  * Rechtsman  _et al._ [2006] M. Rechtsman, F. Stillinger, and S. Torquato, Designed interaction potentials via inverse methods for self-assembly, Physical Review E 73, 011406 (2006).
  * Ma and Ferguson [2019] Y. Ma and A. L. Ferguson, Inverse design of self-assembling colloidal crystals with omnidirectional photonic bandgaps, Soft matter 15, 8808 (2019).
  * King  _et al._ [ress] E. M. King, C. X. Du, Q.-Z. Zhu, S. S. Schoenholz, and M. P. Brenner, Programmable patchy particles for materials design, Proceedings of the National Academy of Sciences (in press).
  * Bradbury  _et al._ [2018] J. Bradbury, R. Frostig, P. Hawkins, M. J. Johnson, C. Leary, D. Maclaurin, G. Necula, A. Paszke, J. VanderPlas, S. Wanderman-Milne, and Q. Zhang, JAX: composable transformations of Python+NumPy programs (2018).
  * Jung  _et al._ [2023] Y.-j. Jung, H. Kim, H.-K. Cheong, and Y.-b. Lim, Magnetic control of self-assembly and disassembly in organic materials, Nature Communications 14, 3081 (2023).
  * Tottori  _et al._ [2013] S. Tottori, L. Zhang, K. E. Peyer, and B. J. Nelson, Assembly, disassembly, and anomalous propulsion of microscopic helices, Nano letters 13, 4263 (2013).
  * Martínez-Pedrero  _et al._ [2020] F. Martínez-Pedrero, F. Ortega, J. Codina, C. Calero, and R. G. Rubio, Controlled disassembly of colloidal aggregates confined at fluid interfaces using magnetic dipolar interactions, Journal of colloid and interface science 560, 388 (2020).
  * Kostiainen  _et al._ [2011] M. A. Kostiainen, P. Ceci, M. Fornara, P. Hiekkataipale, O. Kasyutich, R. J. Nolte, J. J. Cornelissen, R. D. Desautels, and J. van Lierop, Hierarchical self-assembly and optical disassembly for controlled switching of magnetoferritin nanoparticle magnetism, Acs Nano 5, 6394 (2011).
  * Zhang and Glotzer [2004] Z. Zhang and S. C. Glotzer, Self-assembly of patchy particles, Nano letters 4, 1407 (2004).
  * Glotzer and Solomon [2007] S. C. Glotzer and M. J. Solomon, Anisotropy of building blocks and their assembly into complex structures, Nature materials 6, 557 (2007).
  * Rogers [2020] W. B. Rogers, A mean-field model of linker-mediated colloidal interactions, The Journal of chemical physics 153 (2020).
  * Rogers and Crocker [2011] W. B. Rogers and J. C. Crocker, Direct measurements of dna-mediated colloidal interactions and their quantitative modeling, Proceedings of the National Academy of Sciences 108, 15687 (2011).
  * Greener [2023] J. G. Greener, Differentiable simulation to develop molecular dynamics force fields for disordered proteins, bioRxiv , 2023 (2023).
  * Engel  _et al._ [2022] M. C. Engel, J. A. Smith, and M. P. Brenner, Optimal control of nonequilibrium systems through automatic differentiation, arXiv preprint arXiv:2201.00098 (2022).
  * Sanchez-Gonzalez  _et al._ [2020] A. Sanchez-Gonzalez, J. Godwin, T. Pfaff, R. Ying, J. Leskovec, and P. Battaglia, Learning to simulate complex physics with graph networks, in  _International conference on machine learning_ (PMLR, 2020) pp. 8459–8468.
  * Allen  _et al._ [2022] K. R. Allen, T. Lopez-Guevara, K. Stachenfeld, A. Sanchez-Gonzalez, P. Battaglia, J. Hamrick, and T. Pfaff, Physical design using differentiable learned simulators, arXiv preprint arXiv:2202.00728 (2022).
  * Bar-Sinai  _et al._ [2019] Y. Bar-Sinai, S. Hoyer, J. Hickey, and M. P. Brenner, Learning data-driven discretizations for partial differential equations, Proceedings of the National Academy of Sciences 116, 15344 (2019).
  * Baydin  _et al._ [2018] A. G. Baydin, B. A. Pearlmutter, A. A. Radul, and J. M. Siskind, Automatic differentiation in machine learning: a survey, Journal of Marchine Learning Research 18, 1 (2018).