# Nonparametric Identification and Estimation of Earnings Dynamics using a
Hidden Markov Model: Evidence from the PSID

Tong Zhou  Department of Computer Science
Johns Hopkins University
Baltimore, United States
Email<EMAIL_ADDRESS>
DOI: 10.1109/ICAIBD57115.2023.10206080

###### Abstract

This paper presents a hidden Markov model designed to investigate the complex
nature of earnings persistence. The proposed model assumes that the residuals
of log-earnings consist of a persistent component and a transitory component,
both following general Markov processes. Nonparametric identification is
achieved through spectral decomposition of linear operators, and a modified
stochastic EM algorithm is introduced for model estimation. Applying the
framework to the Panel Study of Income Dynamics (PSID) dataset, we find that
the earnings process displays nonlinear persistence, conditional skewness, and
conditional kurtosis. Additionally, the transitory component is found to
possess non-Gaussian properties, resulting in a significantly asymmetric
distributional impact when high-earning households face negative shocks or
low-earning households encounter positive shocks. Our empirical findings also
reveal the presence of ARCH effects in earnings at horizons ranging from 2 to
8 years, further highlighting the complex dynamics of earnings persistence.

###### Index Terms:

Hidden Markov Model, Panel Data, Nonparametric Identification, Modified
Stochastic EM, PSID

##  I Introduction

Earnings dynamics is a fascinating and important area in economics, with
significant implications for understanding economic agents’ consumption
decisions. Macroeconomists employ life-cycle models and profiles of agents’
earnings dynamics to examine their various responses within the economy,
laying the foundation for the creation of sensible policies to manage business
cycles. In a broader context, the nature of earnings dynamics is crucial in
addressing a wide range of economic issues, including income inequality,
optimal design of fiscal policies and insurance programs, economic mobility,
and human capital development. As such, accurately characterizing earnings
dynamics enables more effective management and a deeper understanding of a
country’s economy.

We utilize a parsimonious specification of the earnings process, where log-
earnings consist of an unobserved persistent shock and an unobserved
transitory shock. The literature on earnings process specifications varies in
its focus on the distinction between these two types of shocks, a concept that
can be traced back to Nobel laureate Milton Friedman’s renowned permanent
income hypothesis (PIH). Although there are numerous models of earnings
dynamics, most tend to concentrate on linear specifications for these two
hidden components, inherently excluding the possibility of nonlinear
transmission of earnings shocks.

In this paper, we introduce a new nonparametric framework to explore earnings
dynamics. Both the persistent and transitory components are modeled as two
generic first-order Markov processes. Apart from the first-order restriction,
no further assumptions are imposed on the model. In essence, our specification
establishes a hidden Markov model (HMM) with two latent state variables. Our
focus is on identifying the two Markov kernels, specifically, the conditional
distributions of the persistent component given its past and the conditional
distribution of the transitory component given its past.

We propose a two-step stochastic EM algorithm for estimating the model. In the
E-step, we use an MCMC procedure to obtain draws for the two hidden components
through a likelihood-based approach. In the M-step, we perform a maximization
procedure on a series of quantile regressions with imputed values for hidden
covariates. The iteration continues until the expected likelihood is
maximized.

##  II Materials and Methods

###  II-A Model

####  II-A1 Setup

In line with the conventions of earnings dynamic literature, we use $\log Y$
to represent the real (log) earnings, and it can be decomposed into the
explanatory part, a persistent component $U$ and a transitory component $V$.
The earnings process for each household $i$ at time $t$ is as follows:

$\log(Y_{it})=\mathbf{z}_{it}^{\prime}\bm{\beta}+U_{it}+V_{it},$ (1)

where $\mathbf{z}_{it}$ is a set of observed demographics and known by agents
at $t$. We let $y_{it}=\log(Y_{it})-\mathbf{z}_{it}^{\prime}\bm{\beta}$ denote
the log of real income net of predictable individual components.

We assume both $U_{it}$ and $V_{it}$ follow some general unknown functions
$H_{t}(U_{i,t-1},\eta_{it})$ and $Q_{t}(V_{i,t-1},\varepsilon_{it})$, where
$\eta_{it}$ and $\varepsilon_{it}$ are assumed to follow conditional standard
uniform distributions, i.e.

$\displaystyle\eta_{it}|(U_{i,t-1},U_{i,t-2},\cdots)$
$\displaystyle=\mathsf{Unif}(0,1),~{}~{}t=2,\cdots,T$ (2)
$\displaystyle\varepsilon_{it}|(V_{i,t-1},V_{i,t-2},\cdots)$
$\displaystyle=\mathsf{Unif}(0,1),~{}~{}t=2,\cdots,T.$ (3)

This general nonparametric setting offers greater flexibility for studying the
persistence of earnings dynamics and encompasses many earnings dynamic models
as special cases including the canonical earnings dynamics models, where the
persistent component follows a unit-root process. Given that both processes
are unobserved, a Bernoulli instrumental variable $\omega(C_{it})$ is required
to differentiate them, where $\omega(\cdot)$ is a known transformation of
consumption data $C_{it}$ for agent $i$ at $t$. Since the purpose of
$\omega(C_{it})$ is merely to distinguish the two Markov kernels, it suffices
for our purposes to use a logistic function, i.e.
$\mathbb{P}(\omega(C_{it})=1|U_{it})=1/(1+\exp(-\beta_{0}-\beta_{1}U_{it}))$.
The rationale behind this setup can be found in the works of [3, 4, 5].

Putting above discussions together, we have the complete model setup

$\displaystyle y_{it}$ $\displaystyle=$ $\displaystyle U_{it}+V_{it}$
$\displaystyle U_{it}$ $\displaystyle=$ $\displaystyle
H_{t}(U_{i,t-1},\eta_{it})$ $\displaystyle V_{it}$ $\displaystyle=$
$\displaystyle Q_{t}(V_{i,t-1},\varepsilon_{it})$
$\displaystyle\mathbb{P}(\omega(C_{it})=1)$ $\displaystyle=$ $\displaystyle
1/(1+\exp(-\beta_{0}-\beta_{1}U_{it}))$

####  II-A2 Assumptions

We will outline the assumptions needed to identify the model. Our
identification strategy relies on the powerful spectral decomposition of
linear operators. A thorough overview of this approach can be found in the
work of [11].

######  Assumption 1.

  1. 1.

(First-order Markov) Both $U_{it}$ and $V_{it}$ follow a generic first-order
Markov process;

  2. 2.

(Conditional uniform distribution) Both $\eta_{it}$ and $\varepsilon_{it}$
follow conditional standard uniform distributions

  3. 3.

(Monotonicity) The unknown condition quantile function $\tau\mapsto
H_{t}(U_{i,t-1},\tau)$ and $\tau\mapsto Q_{t}(V_{i,t-1},\tau)$ are strictly
increasing for $\tau\in(0,1)$.

  4. 4.

(Invertibility) The conditional distribution functions $F(U_{it}|U_{i,t-1})$
and $F(V_{it}|V_{i,t-1})$ are both invertible w.r.t. their respective
arguments $U_{it}$ and $V_{it}$ for each $i$ and $t$.

Assumption 1) states that $U_{it}$ and $V_{it}$ have only one-period memory of
their past. This condition imposes dynamic exclusion restrictions that aid in
obtaining nonparametric identifications. This assumption is also commonly made
in structural economic models. Although it can be relaxed to allow for higher-
order Markov process, we maintain the first-order Markovian assumption in this
paper for simplicity. Assumption 2) normalizes the error terms $\eta_{it}$ and
$\varepsilon_{it}$ to follow standard uniform distributions. This setup
enables us to discuss consequences of shocks along the rank of $U_{i,t-1}$ and
$V_{i,t-1}$. This representation also nests the canonical model of earnings
dynamics as a special case where $U_{it}$ is assumed to follow a unit-root
process, i.e.

$U_{i,t+1}=U_{it}+\nu_{i,t+1}$

where $\nu_{i,t+1}=F^{-1}(\eta_{i,t+1})$ is the inverse function of the CDF of
$\eta_{i,t+1}$. Assumption 3) guarantees that $U_{it}$ and $V_{it}$ have
absolutely continuous distributions. Assumption 1) - 3) combined imply that
for all $\tau\in(0,1)$, $H_{t}(U_{i,t-1},\tau)$ happens to be the
$\tau$-conditional quantile of $U_{it}$ given $U_{i,t-1}$. This relationship
also holds for $Q_{t}(V_{i,t-1},\tau)$. Assumption 4) is furnished to
facilitate identification of the nonlinear functions $H_{t}$ and $Q_{t}$. The
monotonicity restriction on $H_{t}$ and $Q_{t}$ are necessary for the
existence of their marginal densities $f(U_{it}|U_{i,t-1})$ and
$f(V_{it}|V_{i,t-1})$. It is not a sufficient condition because a stronger
condition of absolute continuity on the distribution function
$F_{V_{t}|V_{t-1}}$ cannot be weakened. However, since it is rare that a
distribution function is continuous but not absolutely continuous, assumption
4) can be almost equivalent to the existence of the two marginal densities.

######  Assumption 2 (independence).

Two random vectors $(\eta_{i2},\cdots,\eta_{iT},U_{i1})$ and
$(\varepsilon_{i2},\cdots,\varepsilon_{iT},V_{i1})$ are statistically
independent. The Bernoulli random variable $\omega(C_{t})$ is independent of
$V_{t}$ for all $t$.

This assumption suggests that the persistent process
$\left\\{U_{it}\right\\}_{t=1}^{T}$ and the transitory process
$\left\\{V_{it}\right\\}_{t=1}^{T}$ are statistically independent. This
restriction allows for the common deconvolution technique of separating two
unknown probability densities. For instance, once one of the marginal
densities $f(U_{it})$ or $f(V_{it})$ is identified, the other one will also be
automatically identified through the deconvolution argument.

Since our identification strategy relies on the technique of manipulating
linear operators, we provide the definition of linear operator here to
facilitate our later discussions. Let $\mathcal{L}^{p}(F_{U})$ denote the
collection of functions of variable $U$ for which its $p$-th moment is finite,
i.e. $g\in\mathcal{L}^{p}(F_{u})$ implies

$\|g\|_{\mathcal{L}^{p}}=\left(\int_{\mathcal{U}}g(u)\mathrm{d}F_{U}(u)\right)^{\frac{1}{p}}<\infty,$
(4)

where $\mathcal{U}$ denotes the support of $U$. The definition for the space
$\mathcal{L}^{q}(F_{V})$ is similar.

Now we define a linear operator

$\mathcal{L}_{V|U}:\mathcal{L}^{p}(F_{U})\to\mathcal{L}^{q}(F_{V}),$ (5)

where $p,q\geq 1$. Specifically, for any $g\in\mathcal{L}^{p}(F_{U})$, we have

$\mathcal{L}_{V|U}g=\int_{\mathcal{U}}f_{V|U}(v|u)g(u)\mathrm{d}u\in\mathcal{L}^{q}(F_{V}),$
(6)

where the function $f_{V|U}$ is called the kernel of the linear operator
$\mathcal{L}_{V|U}$. This expression is particularly useful when multiple
linear operators are present, since we do not need to introduce new notations
for each involved operator.

######  Assumption 3.

There exist variables $Y_{it}$ such that

  1. 1.

For any $y_{t}$ can $\widetilde{c}_{t}$, there exists a $y_{t-1}$ and
$\widetilde{c}_{t-2}$ and a neighborhood $\mathcal{N}^{r}$ around
$(y_{t},\widetilde{c}_{t-1},y_{t-1},\widetilde{c}_{t-2})$ such that, for any
$(y_{t}^{\prime},\tilde{c}_{t-1}^{\prime},y_{t-1}^{\prime},\tilde{c}_{t-2}^{\prime})\in\mathcal{N}^{r}$,
the linear operator
$\mathcal{L}_{Y_{t-2},y_{t-1}^{\prime},\tilde{c}_{t-2}^{\prime},y_{t}^{\prime},\tilde{c}_{t-1}^{\prime},Y_{t+1}}$
is one-to-one.

  2. 2.

For nay $y_{t}$ and $\tilde{c}_{t-1}$, the linear operator
$\mathcal{L}_{Y_{t=1}|y_{t},\tilde{c}_{t-1},U_{t-1},V_{t}}$ is one-to-one.

  3. 3.

For any $y_{t-1}$ and $\tilde{c}_{t-2}$, the linear operator
$\mathcal{L}_{Y_{t-2},y_{t-1},\tilde{c}_{t-2},Y_{t}}$ is one-to-one.

######  Assumption 4.

  1. 1.

The characteristic function of $(U_{i1},\cdots,U_{iT})$ and
$(V_{i1},\cdots,V_{iT})$ do not vanish on the real line.

  2. 2.

The characteristic function of $(U_{i1},\cdots,U_{iT})$ and
$(V_{i1},\cdots,V_{iT})$ are absolutely continuous.

Assumption 4.1) is commonly made to achieve nonparametric identification (see
[6]). For univariate distributions, this assumption rules out certain families
of distributions, e.g., truncated normal, symmetric uniform and many discrete
distributions. Assumption 4.2) is made to facilitate the deconvolution
argument and also implies that the joint distributions of
$(U_{i1},\cdots,U_{iT})$ and $(V_{i1},\cdots,V_{iT})$ exist.

To avoid cluttered notations, we simplify the notations by omitting the
subscript $i$ without causing confusions. In the following derivations, we
define $\widetilde{C}_{t-1}:=\omega(C_{t-1})$.

######  Assumption 5 (Uniqueness of spectral decomposition).

For any $(Y_{t},\widetilde{C}_{t-1})$ and any
$(u_{t-1},v_{t})\neq(u_{t-1}^{\prime},v_{t}^{\prime})$, there exists a
$(y_{t-1},\widetilde{c}_{t-2})$ and corresponding neighborhood
$\mathcal{N}^{r}$ satisfying Assumption 3.1), such that for some
$(y_{t}^{\prime},\tilde{c}_{t-1}^{\prime},y_{t-1}^{\prime},\tilde{c}_{t-2})\in\mathcal{N}^{r}$
with $(y_{t}^{\prime},\tilde{c}_{t-1})\neq(y_{t},\tilde{c}_{t-1})$ and
$(y_{t-1}^{\prime},\tilde{c}_{t-2})\neq(y_{t-1},\tilde{c}_{t-2})$:

$0<k(y_{t},\tilde{c}_{t-1},y_{t-1}^{\prime},\tilde{c}_{t-2},y_{t-1},\tilde{c}_{t-2},u_{t-1},v_{t})<C<\infty$

and

$\displaystyle
k(y_{t},\tilde{c}_{t-1},y_{t-1}^{\prime},\tilde{c}_{t-2},y_{t-1},\tilde{c}_{t-2},u_{t-1},v_{t})\neq$
$\displaystyle
k(y_{t},\tilde{c}_{t-1},y_{t-1}^{\prime},\tilde{c}_{t-2},y_{t-1},\tilde{c}_{t-2},u_{t-1}^{\prime},v_{t}^{\prime})$

where

$\displaystyle
k(y_{t},\tilde{c}_{t-1},y_{t-1}^{\prime},\tilde{c}_{t-2},y_{t-1},\tilde{c}_{t-2},u_{t-1},v_{t})=$
$\displaystyle\frac{f(y_{t},\tilde{c}_{t-1}|y_{t-1},\tilde{c}_{t-2},u_{t-1},v_{t})f(y_{t}^{\prime},\tilde{c}_{t-1}^{\prime}|y_{t-1}^{\prime},\tilde{c}_{t-2}^{\prime},u_{t-1},v_{t}))}{f(y_{t}^{\prime},\tilde{c}_{t-1}^{\prime}|y_{t-1},\tilde{c}_{t-2},u_{t-1},v_{t}f(u_{t-1},v_{t}|y_{t-1}^{\prime},\tilde{c}_{t-2}^{\prime},u_{t-1},v_{t}))}$

######  Assumption 6 (normalization).

The Markov kernels are normalized by $\mathbb{E}[U_{t+1}|U_{t}]=U_{t}$ and
$\mathbb{E}[V_{t+1}|V_{t}]=0$.

In the eigenfunctions $f(y_{t+1}|y_{t},\tilde{c}_{t-1},u_{t-1},v_{t})$, both
$u_{t-1}$ and $v_{t}$ are unobserved and continuously distributed. Assumption
6 is made to differentiate and identify the two components.

######  Assumption 7 (Stationarity).

For any $2\leq t\leq T$, the Markov kernels is time-invariant, i.e.,

$\displaystyle
f(Y_{t},\widetilde{C}_{t-1},U_{t-1},V_{t}|Y_{t-1},\widetilde{C}_{t-2},U_{t-2},V_{t-1})$
$\displaystyle=f(Y_{3},\widetilde{C}_{2},U_{2},V_{3}|Y_{2},\widetilde{C}_{1},U_{1},V_{2})$

Assumption 7 is not necessary for identification of the Markov density. It
eases our derivations. From the next section, we can see that only five
periods of data is sufficient for achieving nonparametric identification.

###  II-B Identification

Based on the assumptions made in the previous section, identification can be
accomplished by applying Theorem 9 from [11] in dynamic settings. This
strategy can be better understood in Figure 1, where the dependence structures
and dynamic exclusion restrictions can be easily visualized.

$(Y_{t-2},\widetilde{C}_{t-3})$$(Y_{t-1},\widetilde{C}_{t-2})$$(Y_{t},\widetilde{C}_{t-1})$$(Y_{t+1},\widetilde{C}_{t})$$U_{t-2}$$U_{t-1}$$~{}U_{t~{}~{}}$$U_{t+1}$$V_{t-2}$$V_{t-1}$$~{}V_{t~{}~{}}$$V_{t+1}$
Figure 1: Graphical illustration of earnings dynamics

We state the main identification theorem

######  Theorem 1 (Identification).

Under Assumption 1 - Assumption 7, the density

$f(Y_{t+1},\widetilde{C}_{t},Y_{t},\widetilde{C}_{t-1},Y_{t-1},\widetilde{C}_{t-2},Y_{t-2},\widetilde{C}_{t-3})$

for any $t\in\left\\{4,\cdots,T-1\right\\}$ uniquely determines the densities
$f(Y_{t},\widetilde{C}_{t-1},U_{t-1},V_{t}|Y_{t-1},\widetilde{C}_{t-2},U_{t-2},V_{t-1})$.

Theorem 1 implies that our interests of Markov kernels can be identified by
basic probability rules, the Bayes rules and the deconvolution technique.

######  Corollary 1.

Under Assumption 1 - Assumption 7, the Markov kernels $f_{V_{t}|V_{t-1}}$,
$f_{U_{t}|U_{t-1}}$ and marginal distributions $f_{U_{t}}$ and $f_{V_{t}}$ are
uniquely identified, for $t=4,\dots,T-1$.

###  II-C Estimation

We introduce a modified stochastic EM algorithm (MSEM) to estimate this HMM,
while the stochastic EM was originally proposed by [10]. The MSEM provides a
much more faster implementation of the estimation by replacing the likelihood
with the objective functions of quantile regression models. The MSEM is
similar to the one presented in [1]. The difference lies in the fact that
their paper involves only one state variables. Specifically, for any
$\tau\in(0,1)$ we employ the following estimating equations

$\displaystyle U_{it}$
$\displaystyle=\sum_{k=0}^{K_{1}}a_{k}^{H}(\tau)\phi_{k}(U_{i,t-1},\mathrm{age}_{it})$
$\displaystyle V_{it}$
$\displaystyle=\sum_{k=0}^{K_{2}}a_{k}^{Q}(\tau)\phi_{k}(V_{i,t-1},\mathrm{age}_{it})$
$\displaystyle U_{i1}$
$\displaystyle=\sum_{k=0}^{K_{3}}a_{k}^{H_{1}}(\tau)\phi_{k}(\mathrm{age}_{i1})$
$\displaystyle V_{i1}$
$\displaystyle=\sum_{k=0}^{K_{2}i}a_{k}^{Q_{1}}(\tau)\phi_{k}(\mathrm{age}_{i1}),$

where $\phi_{k}$ is the Hermite polynomials. We selected different orders of
polynomials for the four equations to maximize the likelihood. The quantile-
based estimation strategy provides a flexible specification of the Markov
kernels. [2] applies this estimation strategy to estimate the smoking effects
of women during pregnancy on children’s birthweights.

Another advantage of using quantile-based estimation is that the original
nonparametric estimation problem is reduced to estimating a finite number of
parameters, i.e., the coefficients of the Hermite polynomials. We discuss
$U_{it}$ as an example: the functions $a_{k}^{H}(\tau)$ are modeled as
piecewise-polynomial interpolating splines on equi-length intervals
$[\tau_{1},\tau_{2}],[\tau_{2},\tau_{3}],\cdots,[\tau_{I-1},\tau_{I}]$ that
partition the unit interval $(0,1)$. In other words, we need to estimate
$a_{k}^{H}(\tau)$ for each interval of $\tau$ and $k$. Additionally, the
objective function of quantile regressions can be used as a surrogate
likelihood. Since it is a convex function, the implementation can be fast.
Once those $a_{\tau}^{H}$ ’s are obtained, we are finished with the estimation
of $U_{it}$.

We still take $U_{it}$ as an example to illustrate the MSEM algorithm. We
start with an initial value for the parameter vector $\widehat{\theta}^{(0)}$.
Each iteration follows the following two steps until convergence of the
$\widehat{\theta}^{(s)}$ in the $s$-th iteration:

  * •

_Stochastic E-step:_ Draw $U_{i}^{(m)}=(U_{i1}^{(m)},\cdots,U_{iT}^{(m)})$ for
$m=1,\cdots,M$ from $f_{i}(\cdot;\widehat{\theta}^{(s)})$.

  * •

_M-step:_ Compute

$\widehat{\theta}^{(s+1)}=\operatorname*{arg\,min}_{\theta}\sum_{i=1}^{N}\sum_{m=1}^{M}R(y_{i},U_{i}^{(m)};\theta),$

where $R(\cdot)$ is the surrogate likelihood, i.e. the objective function of
the piece-wise quantile regressions. In the E-step, we use a random-walk
Metropolis-Hastings algorithm for drawing $U_{i}^{(m)}$ in the E-step. The
M-step consists of a number of quantile regressions. For instance, for each
$\ell$, the parameters $a_{k}^{H}(\tau_{\ell})$ are updated as

$\displaystyle\min_{(a_{0\ell}^{H},\cdots,a_{K\ell}^{H})}\sum_{i=1}^{N}$
$\displaystyle\sum_{t=2}^{T}\sum_{m=1}^{M}$
$\displaystyle\rho_{\tau_{\ell}}\left(U_{it}^{(m)}-\sum_{k=0}^{K}a_{k\ell}^{H}\varphi_{k}\left(U_{i,t-1}^{(m)},\textsf{age}_{it}\right)\right),$

where $\rho_{\tau}(u)=u(\tau-\bm{1}(u\leq 0))$ is the check function in
standard quantile regressions and $\bm{1}(\cdot)$ denotes the indicator
function, first introduced by [8].

[10] examined the statistical properties of the stochastic EM algorithm within
a likelihood case. He provided certain conditions under which the Markov chain
$\widehat{\theta}^{(s)}$ is ergodic. He also outlined the asymptotic
distribution of $\widehat{\theta}$. [1] characterized the asymptotic
distribution of $\widehat{\theta}$ in a manner that aligns with our model,
specifically when utilizing the surrogate likelihood during the M-step.

The MSEM algorithm can be summarized as follows

##  III Results and Discussions

###  III-A Data

As the longest running household panel survey in the world, the PSID dataset
contains a large amount of data in the US over 50 years. The dataset includes
a wide range of variables on income, employment, education, economic, social,
health-related factors and many other aspects of life for each individual and
their family members. However, the size, variety and complexity of the dataset
make it challenging to analyze using traditional statistical methods.

The second challenge comes from the longitudinal nature, which means that it
follows the same individuals and families over time. This presents unique
challenges for analysis, such as handling missing data, attrition, and changes
in the variables of interest over time.

The third one is data quality. The quality of PSID can vary over time, as
changes in survey methodology or sample composition can affect the accuracy
and reliability of the data. This requires careful attention to data cleaning
and quality control procedures.

Our study draws on the PSID as our primary data source. To ensure high quality
data and better compare our empirical findings with other literature, our
sample selection procedure and data preprocessing mainly follow the works of
[3, 4].

###  III-B Empirical Findings

####  III-B1 Densities and Moments

Figure 2 illustrates the marginal distributions of the persistent and
transitory earnings components at the mean age. The persistent component
$U_{it}$ displays small deviations from Gaussianity. However, the marginal
distribution of $V_{it}$ provides strong evidence to reject Gaussianity owing
to its high kurtosis and fat tails. It is worth noting that the density of
$V_{it}$ in our model is less spiky than that in [1]. A possible explanation
for this difference could be the mutual dependence structure of $V_{it}$,
which is governed by its first-order Markovian property, whereas they are
assumed to be mutually independent across $t$ in their paper. In Figure 3, we
report the conditional sknewness for $\tau=11/12$, for both the $U$ component
and the $V$ component. These two panels show a similar pattern:
$U_{it}$($V_{it}$) is positively skewed for low values of
$U_{i,t-1}$($V_{i,t-1}$), and negatively skewed for high values of
$U_{i,t-1}$($V_{i,t-1}$).

Figure 2: Marginal distributions of persistent and transitory earnings
components

Figure 3: Conditional skewness of $U$ component and $V$ component.

####  III-B2 Nonlinear persistence

This experiment examines the marginal effects of persistent and transitory
shocks. For linear specifications of earnings dynamics, their marginal effects
would be constant by design. In the canonical model of earnings dynamics, for
example, where the innovation is a random walk, then the marginal effect is
$1$, regardless of $U_{i,t-1}$ and $\tau$. In contrast, our model allows the
persistence of $U_{i,t-1}$ to depend on the magnitude and direction of the
shock. As a result, the persistence of a shock to $U_{i,t-1}$ depends on the
size and sign of current and future shocks. In particular, our model enables
specific shocks to erase the memory of past shocks. Furthermore, the
interaction between the shock $\eta_{it}$ and the lagged persistent component
$U_{i,t-1}$ is a central feature of our nonlinear approach. We then estimate
the earnings model, and, given the estimated parameters, we simulate the
model. Figure 5 shows that our nonlinear model reproduces the patterns of
nonlinear persistence well.

Figure 4 indicates the presence of nonlinear persistence, which depends on
both the percentile of past earnings $(\tau_{\textsf{init}})$ and the
percentile of the quantile innovation $(\tau_{\textsf{shock}})$. Figure 6 then
shows the estimated persistence of the earnings component $U_{it}$.
Specifically, the graph shows the marginal effects, evaluated at percentiles
$\tau_{\textsf{init}}$ and $\tau_{\textsf{shock}}$ and at the mean age in the
sample. Persistence in $U$’s is higher than persistence in log-earnings
residuals, consistently with the fact that Figure 6 is net of transitory
shocks. One observation that sets this study apart from [1] is that the
persistence in Figure 6 is higher than 1. For high-earnings households hit by
good shocks and low-earnings households hit by bad shocks, persistence is even
above 1.5, with the persistence in the latter being higher than that in the
former.

Figures 7 and 8 demonstrate that the persistence in $V$’s is generally lower
in magnitude than that in $U$. One notable feature is that when high-earnings
households are hit by bad shocks and low-earnings households are hit by good
shock, the persistence can be negative. The high degree of nonlinearity
displayed here strongly rejects that $V_{t}$ follows an independent process
across $t$. If it did, its nonlinear persistence measures would remain
constant for any $t$.

Figure 4: Estimates of the average derivative of the conditional quantile
function of log-earnings residuals $y_{it}$ given $y_{i,t-1}$ with respect to
$y_{i,t-1}$ in the PSID Figure 5: Estimates of the average derivative of the
conditional quantile function of simulated model. Figure 6: Estimates of the
average derivative of the conditional quantile function of the persistent
component $U$. Figure 7: Estimates of the average derivative of the
conditional quantile function of the transitory component $V$. Figure 8:
Estimates of the average derivative of the conditional quantile function of
the persistent component $V$.

####  III-B3 ARCH effects

Figure 9 presents estimates of log-earnings residuals growth at various
horizons, from 2 to 8 years. All of them suggest the presence of ARCH effects,
which is consistent with findings in the existing literature, such as [9]. The
data also reveal that log-earnings growth is non-Gaussian and displays
negative skewness and high kurtosis. [7] finds similar features in U.S.
administrative data. [1] further highlights the sknewness and excess kurtosis
of log-earnings growth at long horizons are primarily due to the non-
Gaussianity of the transitory component.

Figure 9: Densities of log-earnings growth at various horizons.

##  IV Conclusion

We develop a nonparametric identification strategy for modeling earnings
dynamics, differentiating the two unobserved components based on their
distinct impact on household consumption. We also propose an modified
stochastic EM algorithm for estimating this model. The identification tool
relies on the assumptions that several linear operators are one-to-one.

In analyzing PSID, the empirical results reveal notable nonlinearities in both
persistence component and transitory component. Specifically, substantial
nonlinear persistence and conditional skewness are observed in both
components. These findings suggest that the earnings shocks to a household
depend on both the history of past shocks and the household’s past relative
wealth. In particular, persistence is higher for high-earnings households hit
by good shocks and low-earnings households hit by bad shocks, while it is
lower for high-earnings household hit by bad shocks and low-earnings
households hit by good shocks. These features align with similar observations
in the PSID that previous earnings dynamic models cannot capture. We also find
some other features such as ARCH effects that have been documented in other
literature.

## Acknowledgment

The authors would like to thank the University of Michigan for providing the
Panel Study of Income Dynamics (PSID) data, which was essential to the success
of our research. We are also deeply grateful to Professor Roger Koenker for
his valuable suggestions and clarifications on issues about quantile
regression models.

## References

  * [1] Arellano, Manuel, Richard Blundell, and Stéphane Bonhomme. “Earnings and consumption dynamics: a nonlinear panel data framework.” Econometrica 85.3 2017: 693-734.
  * [2] Arellano, Manuel, and Stéphane Bonhomme. “Nonlinear panel data estimation via quantile regressions.” Econometric Theory. 2016: C61-C94.
  * [3] Blundell, Richard, Luigi Pistaferri, and Ian Preston. “Consumption inequality and partial insurance.” American Economic Review 98.5 2008: 1887-1921.
  * [4] Blundell, Richard, Luigi Pistaferri, and Itay Saporta-Eksten. “Consumption inequality and family labor supply.” American Economic Review 106.2 2016: 387-435.
  * [5] Blundell, Richard, Luigi Pistaferri, and Itay Saporta-Eksten. “Children, time allocation, and consumption insurance.” Journal of Political Economy 126.S1 2018: S73-S115.
  * [6] D’Haultfoeuille, Xavier. “On the completeness condition in nonparametric instrumental problems.” Econometric Theory 27.3 2011: 460-471.
  * [7] Guvenen F, Karahan F, Ozkan S, Song J. “What do data on millions of US workers reveal about life-cycle earnings risk?”. National Bureau of Economic Research; 2015 Feb 2.
  * [8] Koenker R, Bassett Jr G. “Regression quantiles”. Econometrica: journal of the Econometric Society. 1978 Jan 1:33-50.
  * [9] Meghir C, Pistaferri L. “Income variance dynamics and heterogeneity”. Econometrica. 2004 Jan;72(1):1-32.
  * [10] Nielsen SF. The stochastic EM algorithm: estimation and asymptotic results. Bernoulli. 2000 Jun 1:457-89.
  * [11] Schennach S. “Measurement systems. Journal of Economic Literature”. 2022 Dec;60(4):1223-63.


#  Physically Compatible 3D Object
Modeling from a Single Image

Minghao Guo MIT CSAIL  Bohan Wang Pingchuan Ma MIT CSAIL  Tianyuan Zhang MIT
CSAIL  Crystal Elaine Owens MIT CSAIL  Chuang Gan UMass Amherst  MIT-IBM
Waston AI Lab  Joshua B. Tenenbaum MIT CSAIL  MIT BCS  Kaiming He MIT CSAIL
Wojciech Matusik MIT CSAIL

###### Abstract

We present a computational framework that transforms single images into 3D
physical objects. The visual geometry of a physical object in an image is
determined by three orthogonal attributes: mechanical properties, external
forces, and rest-shape geometry. Existing single-view 3D reconstruction
methods often overlook this underlying composition, presuming rigidity or
neglecting external forces. Consequently, the reconstructed objects fail to
withstand real-world physical forces, resulting in instability or undesirable
deformation – diverging from their intended designs as depicted in the image.
Our optimization framework addresses this by embedding physical compatibility
into the reconstruction process. We explicitly decompose the three physical
attributes and link them through static equilibrium, which serves as a hard
constraint, ensuring that the optimized physical shapes exhibit desired
physical behaviors. Evaluations on a dataset collected from Objaverse
demonstrate that our framework consistently enhances the physical realism of
3D models over existing methods. The utility of our framework extends to
practical applications in dynamic simulations and 3D printing, where adherence
to physical compatibility is paramount.

**footnotetext: Corresponding author.

Figure 1:  Existing methods for single-view reconstruction often result in
objects that, when subjected to real-world physical forces (such as gravity)
and user-required mechanical materials, exhibit problematic behaviors such as
toppling over (top left) and undesirable deformation (top right), diverging
from their intended depiction in the input images. In contrast, our approach
produces physical objects that maintain stability (bottom left) and mirror the
objects’ static equilibrium state captured in the input images (bottom right).

##  1 Introduction

The field of single-image 3D shape modeling has experienced significant
advancements over the past years, largely propelled by advances in single-view
reconstruction techniques. These methods, ranging from generating multi-view
consistent images for per-scene 3D reconstruction [22, 23, 21, 24, 30, 20], to
employing large reconstruction models (LRMs) for feedforward inference [12,
39, 42, 48, 45, 38], have enhanced the geometric quality and visual fidelity
of the 3D shapes to unprecedented levels.

However, reconstructing a 3D shape from an image often aims to be beyond a
mere visualization. These generated objects find applications in virtual
environments such as filming and gaming, as well as in tangible fields like
industrial design and engineering. Despite their diverse applications, a
common oversight in many current single-view reconstruction methods is the
negligence of physical principles. As shown in the top row of Fig. 1, when
subjected to real-world physics such as gravity, these 3D objects produced
from these techniques exhibit issues such as instability and undesired
deformation, diverging from their depiction in the input images. Such
inconsistency can significantly undermine the practical utility of the models,
as they fail to meet the functional and aesthetic expectations set by the
original image.

Fundamentally, an image is more than a visual representation of an object: It
captures a physical snapshot of the object in a state of static equilibrium,
under the influence of real-world forces. In this context, the geometry seen
in an image is determined by three orthogonal attributes: _mechanical
properties_ , _external forces_ , and _rest-shape geometry_. As shown in the
inset figure, these attributes collectively model the spectrum of potential
static configurations that a physical object might adopt. Reconstructing such
an object from an image is essentially an ill-posed problem, since multiple
combinations of these attributes can result in identical static geometry.
Current methods, however, often overlook this underlying composition; they
typically assume objects are rigid or neglect the impact of external forces.
The reconstructed objects thus merely replicate the visual geometry without
considering the three physical attributes.

To bridge this gap, we explicitly decompose these attributes for
reconstructing a physical object from a single image. Our framework
holistically takes mechanical properties and external forces as predefined
inputs, reflecting typical user specifications in real-world applications like
3D printing and simulations. The output is the rest-shape geometry tailored to
these inputs. These attributes are integrated through the principles of static
equilibrium physics. This explicit decomposition imposes two stringent
physical constraints in object modeling: static equilibrium is enforced as a
_hard constraint_ , and the physical object must conform to user-specified
material properties and external forces. These resulting physical objects are
stable, robust under real-world physics, and are high-fidelity replicas
inferred from the input images, as shown in the bottom row of Fig. 1.

More specifically, we propose _physical compatibility_ optimization, which is
a physically constrained optimization with rest-shape geometry as the
variable. In this setup, the objective is for the modeled physical object to
exhibit desired behaviors, such as matching the geometry depicted in the input
image under external forces and maintaining stability under gravity. The
constraint is the equation of static equilibrium simulation, ensuring that
during optimization, the physical object remains in the equilibrium state,
with internal forces generated by deformation from the rest shape balancing
the external forces. We parameterize the rest-shape geometry using a plastic
deformation field and solve this hard-constrained optimization problem by
using implicit differentiation with gradient descent.

For evaluation, we introduce five metrics designed to comprehensively assess
the physical compatibility of the modeled 3D objects under simulation. These
metrics include image loss between the rendered image of the modeled physical
object and the input image, stability under gravity, as well as measures from
finite element analysis, such as integrity and structural robustness. Our
framework’s versatility is demonstrated by its integration with five distinct
single-view reconstruction methods, each employing unique geometry
representations. Results on a dataset collected from Objaverse [9], consisting
of $100$ shapes, show that our framework consistently produces 3D objects with
enhanced physical compatibility. Furthermore, we demonstrate the practical
utility of our framework through applications in dynamic simulations and 3D
printing fabrication.

#### Related work.

Our method is mainly related to: 1) single-view 3D reconstruction, where our
work emphasizes the integration of physical modeling into the reconstruction
process; 2) physically based 3D modeling, where we incorporate physical
principles as hard constraints within the reconstruction process; and 3)
fabrication-aware shape design, where our work directly constructs physical
objects from a single input image rather than relying on manual creation of
the initial geometry. A comprehensive discussion of related work is provided
in Appendix G.

##  2 Approach

Figure 2: Overall pipeline. Given predefined mechanical properties and
external forces, our pipeline optimizes the rest-shape geometry to ensure that
the shape, when in a state of static equilibrium, aligns with the target image
and meets stability criteria. We visualize the stress distribution of the
static geometry using a colored heat map, illustrating the spatially varying
deformation of the physical object under static equilibrium.

Our objective is to create 3D objects from a single image that are physically
compatible, ensuring that they align with the input image in the static
equilibrium state while also meeting the stability requirements. Governed by
universal physical principles, the physical behavior of an object is
determined by its mechanical properties, external forces, and rest-shape
geometry. Our framework treats the rest-shape geometry as the optimization
variable, assuming that the mechanical properties and external forces are
predefined as inputs. Fig. 2 illustrates the overall pipeline.

###  2.1 Formulation of Physical Compatibility

In our approach, we treat the entity depicted in the input image as a solid
object. We employ Finite Element Method (FEM) for robust solid simulation. The
object is represented by a volumetric mesh, denoted as
$\mathcal{M}=(\mathbf{X},\mathbf{T})$. Here, $\mathbf{X}\in\mathbb{R}^{3N}$
represents the 3D positions of the vertices, with $N$ denoting the total
number of vertices. $\mathbf{T}\in\mathbb{N}^{Z\times K}$ describes the mesh
connectivity, where $Z$ represents the total number of elements and $K$
indicates the number of vertices per element. The mesh in its _rest-shape
geometry_ , which is the state without any internal or external forces
applied, is represented as
$\mathcal{M}_{\mathrm{rest}}=(\mathbf{X}_{\mathrm{rest}},\mathbf{T})$. The
input image depicts the _static geometry_ , which is the deformed geometry of
the object under static equilibrium111Although our implementation employs
_quasi-static equilibrium_ , we use the term _static equilibrium_ across the
paper for consistency., denoted as
$\mathcal{M}_{\mathrm{static}}=(\mathbf{x}_{\mathrm{static}},\mathbf{T})$. In
accordance with Newton’s laws, $\mathbf{x}_{\mathrm{static}}$ adheres to the
following equation:

$\mathbf{f}_{\mathrm{int}}(\mathbf{x}_{\mathrm{static}},\mathbf{X}_{\mathrm{rest}};\Theta)=\mathbf{f}_{\mathrm{ext}}(\mathbf{x}_{\mathrm{static}}),$
(1)

where
$\mathbf{f}_{\mathrm{int}}(\cdot,\cdot;\Theta):\mathbb{R}^{3N}\times\mathbb{R}^{3N}\rightarrow\mathbb{R}^{3N}$
denotes the internal forces exerted by deformed objects transitioning from
$\mathbf{X}_{\mathrm{rest}}$ to $\mathbf{x}_{\mathrm{static}}$,
$\mathbf{f}_{\mathrm{ext}}(\cdot):\mathbb{R}^{3N}\rightarrow\mathbb{R}^{3N}$
embodies the external interaction forces such as gravity, and $\Theta$
represents the mechanical material properties, such as the stiffness of the
object. Eq. 1 reveals that $\Theta$ (mechanical properties),
$\mathbf{f}_{\mathrm{ext}}$ (external forces), and
$\mathbf{X}_{\mathrm{rest}}$ (the rest-shape geometry) collectively determine
the static geometry $\mathbf{x}_{\mathrm{static}}$.

Given $\Theta$ and $\mathbf{f}_{\mathrm{ext}}(\cdot)$, the goal of physically
compatible modeling is to ensure that the rest-shape geometry
$\mathcal{M}_{\mathrm{rest}}$ conforms to given objectives under static
equilibrium. This is formulated as the following optimization problem:

$\displaystyle\min_{\mathbf{X}_{\mathrm{rest}},\mathbf{x}_{\mathrm{static}}}$
$\displaystyle\quad\mathcal{J}(\mathbf{X}_{\mathrm{rest}},\mathbf{x}_{\mathrm{static}})=\mathcal{L}(\mathbf{x}_{\mathrm{static}})+\mathcal{L}_{\mathrm{reg}}(\mathbf{X}_{\mathrm{rest}})$
$\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\mathbf{f}_{\mathrm{int}}(\mathbf{x}_{\mathrm{static}},\mathbf{X}_{\mathrm{rest}};\Theta)=\mathbf{f}_{\mathrm{ext}}(\mathbf{x}_{\mathrm{static}}).$
(2)

Here, $\mathcal{J}(\mathbf{X}_{\mathrm{rest}},\mathbf{x}_{\mathrm{static}})$
is the objective function, consisting of
$\mathcal{L}(\mathbf{x}_{\mathrm{static}})$, which measures the alignment of
the geometry $\mathbf{x}_{\mathrm{static}}$ with the specified target.
$\mathcal{L}_{\mathrm{reg}}(\mathbf{X}_{\mathrm{rest}})$ regularizes the rest-
shape geometry $\mathbf{X}_{\mathrm{rest}}$, with more details discussed in
Section 2.2.

Within the scope of this work, two tasks for
$\mathcal{L}(\mathbf{x}_{\mathrm{static}})$ are considered: 1)
$\mathbf{x}_{\mathrm{static}}$ replicates the geometry depicted in the input
image; and 2) $\mathbf{x}_{\mathrm{static}}$ maintains stability and
inherently remains upright without toppling. In the first scenario, the loss
function is
$\mathcal{L}(\mathbf{x}_{\mathrm{static}})=\|\mathbf{x}_{\mathrm{static}}-\mathbf{X}_{\mathrm{target}}\|_{2}^{2}$
which measures the point-wise Euclidean distance between the static shape and
the target geometry
$\mathcal{M}_{\mathrm{target}}=(\mathbf{X}_{\mathrm{target}},\mathbf{T})$. In
the second scenario, the loss function is
$\mathcal{L}(\mathbf{x}_{\mathrm{static}})=\|\mathrm{proj}_{z}(\mathcal{C}(\mathbf{x}_{\mathrm{static}}))-\hat{\mathcal{C}}\|$,
where $\mathcal{C}(\cdot)$ computes the center of mass of
$\mathcal{M}_{\mathrm{static}}$, $\mathrm{proj}_{z}(\cdot)$ denotes the
projection of the center onto the $z$-plane in world coordinates, and
$\hat{\mathcal{C}}$ represents the target position for the center of mass to
guarantee stability. Minimization of this function ensures the structural
stability of $\mathcal{M}_{\mathrm{static}}$.

It is crucial to highlight that the variables $\mathbf{X}_{\mathrm{rest}}$ and
$\mathbf{x}_{\mathrm{static}}$ are tightly coupled through a hard constraint
in our problem formulation. This constraint, which ensures that the object
remains static equilibrium, is essential to achieving physical compatibility.
Enforcing this configuration guarantees that the 3D physical object conforms
strictly to external forces such as gravity, thereby ensuring the system
adheres to the inherent physical constraints.

###  2.2 Parameterization of Rest-shape Geometry

To solve the optimization problem defined Eq. 2.1, one might consider a
straightforward approach by directly treating $\mathbf{X}_{\mathrm{rest}}$ as
the optimization variable. However, this brings challenges in maintaining the
physical validity of the rest-shape geometry, i.e., there shall be no
inversions or inside-out elements. This non-inversion requirement is typically
enforced through nonlinear inequality constraints [11, 34], leading to
intractable optimization. Drawing inspiration from natural modeling processes
[40], we propose a parameterization of $\mathbf{X}_{\mathrm{rest}}$ by
treating it as the result of plastic deformation applied to an initial
configuration. A _plastic deformation_ can transform objects without the
volume preservation constraint [1]. Specifically, we denote the initial
configuration of the rest-shape geometry as
$\mathcal{M}_{\mathrm{init}}=(\mathbf{X}_{\mathrm{init}},\mathbf{T})$.
$\mathbf{X}_{\mathrm{rest}}$ is implicitly parameterized by the plastic
deformation field $\mathbf{F}_{\mathbf{p}}$ as

$\mathbf{X}_{\mathrm{rest}}:=\phi(\mathbf{F}_{\mathbf{p}};\mathbf{X}_{\mathrm{init}}),\quad\text{with}\quad\mathbf{f}_{\mathrm{int}}(\mathbf{X}_{\mathrm{rest}},\mathbf{X}_{\mathrm{init}};\Theta)=\mathbf{0}.$
(3)

Intuitively, this equation suggests that $\mathbf{X}_{\mathrm{rest}}$ results
from applying plastic strain field $\mathbf{F}_{\mathbf{p}}$ to
$\mathbf{X}_{\mathrm{init}}$ without any external forces. The plastic strain
field $\mathbf{F}_{\mathbf{p}}$ is the collection of transformations, with
each transformation is an $\mathbb{R}^{3\times 3}$ matrix applied to each
material point. Throughout this paper, we also represent plastic deformation
in its vector form as $\mathbf{F}_{\mathbf{p}}\in\mathbb{R}^{9Z}$, which
corresponds to the flattened vector form of the $\mathbb{R}^{3\times 3}$
transformation collection. For a detailed explanation of the computation of
$\mathbf{X}_{\mathrm{rest}}$ from $\mathbf{F}_{\mathbf{p}}$ and its
integration into the static equilibrium, we refer the reader to Appendix B.

There are several benefits using $\mathbf{F_{p}}$ for parameterizing rest-
shape geometry: It exhibits invariance to translation, which ensures that the
spatial positioning of $\mathbf{X}_{\mathrm{init}}$ does not affect the
deformation outcomes. Moreover, the non-inversion requirement can be
efficiently satisfied by constraining the singular values of
$\mathbf{F}_{\mathbf{p}}$, thereby avoiding the need for complicated
inequality constraints. Appendix B provides a comprehensive analysis of these
advantages.

By substituing Eq. 3, we reformulate the optimization problem Eq. 2.1 as
follows:

$\displaystyle\min_{\mathbf{F_{p}},\mathbf{x}_{\mathrm{static}}}$
$\displaystyle\quad\mathcal{J}(\mathbf{F_{p}},\mathbf{x}_{\mathrm{static}})=\mathcal{L}(\mathbf{x}_{\mathrm{static}})+\mathcal{L}_{\mathrm{reg}}(\mathbf{F_{p}})$
$\displaystyle\mathrm{s.t.}$
$\displaystyle\quad\mathbf{f}_{\mathrm{int}}(\mathbf{x}_{\mathrm{static}},\phi(\mathbf{F}_{\mathbf{p}};\mathbf{X}_{\mathrm{init}});\Theta)=\mathbf{f}_{\mathrm{ext}}(\mathbf{x}_{\mathrm{static}}).$
(4)

Here, the optimization variables are $\mathbf{F_{p}}$, where the initial
geometry configuration $\mathbf{X}_{\mathrm{init}}$ is treated as a constant.
The regularization term $\mathcal{L}_{\mathrm{reg}}(\mathbf{F_{p}})$ is
defined as the smoothness of plastic deformation using bi-harmonic energy [5],
represented as
$\mathcal{L}_{\mathrm{reg}}(\mathbf{F_{p}})=\|\mathbf{L}\mathbf{F_{p}}\|_{2}^{2}$,
where $\mathbf{L}\in\mathbb{R}^{9Z\times 9Z}$ denotes the graph Laplacian
matrix, encapsulating the connectivity of the volumetric mesh elements.

###  2.3 Implicit Differentiation-based Optimization

Solving the optimization problem in Eq. 2.2 is non-trivial due to its
nonlinear objective and the nonlinear hard constraint. A straightforward
approach is incorporating the constraint directly into the objective as an
additional loss term; however, this method may lead to imperfect satisfaction
of the constraint, which undermines the fundamental goal of ensuring physical
compatibility.

We resort to implicit differentiation, a technique used in sensitivity
analysis [6], to compute the gradient of the objective function $\mathcal{J}$
with respect to the variable $\mathbf{F_{p}}$. This approach effectively
reduces the dimensionality of the optimization variables since we only need to
calculate the gradient with respect to $\mathbf{F_{p}}$ and also ensures that
the gradient direction takes into account the hard constraint. Specifically,
the gradient is computed as follows:

$\displaystyle\frac{\partial\mathcal{J}}{\partial\mathbf{F_{p}}}=-\left(\frac{\partial\mathcal{L}}{\partial\mathbf{x}_{\mathrm{static}}}\right)\left[\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{x}_{\mathrm{static}}}\right]^{-1}\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{F_{p}}}+\frac{\partial\mathcal{L}_{\mathrm{reg}}}{\partial\mathbf{F_{p}}},$
(5)

where
$\mathbf{f}_{\mathrm{net}}=\mathbf{f}_{\mathrm{int}}-\mathbf{f}_{\mathrm{ext}}$
represents the net forces. A comprehensive derivation of this gradient formula
is provided in Appendix C. By utilizing this gradient, the optimization can be
solved using standard optimization tools, such as the Adam optimizer [17].
This facilitates the integration of our method into existing single-view
reconstruction pipelines.

###  2.4 Implementation Details

Given an input image, we initially utilize off-the-shelf single-view
reconstruction models to obtain the 3D object’s target geometry, ensuring
alignment with the input image. The output of these reconstruction models
varies depending on the geometric representation used. For instance, methods
employing tetrahedral representations, such as TetSphere [11], yields
volumetric meshes that can be directly used as
$\mathcal{M}_{\mathrm{target}}$. Conversely, methods that output surface
meshes [42] or point clouds [38], which are often non-volumetric and typically
non-manifold, require additional processing steps to be suitable for our
computational pipeline. We use TetWild [15], a robust tetrahedral meshing
algorithm, to convert these unstructured outputs into high-quality tetrahedral
meshes, resulting in volumetric mesh $\mathcal{M}_{\mathrm{target}}$. For
initiating the optimization process, we set
$\mathcal{M}_{\mathrm{init}}=\mathcal{M}_{\mathrm{target}}$, assuming that
$\mathcal{M}_{\mathrm{target}}$ is a reasonably good initial approximation for
the optimization. Note that $\mathcal{M}_{\mathrm{init}}$ is not strictly
confined to $\mathcal{M}_{\mathrm{target}}$; any volumetric mesh could
potentially serve as the initial approximation, given the flexibility of
$\mathbf{F_{p}}$ to accommodate spatially varying deformations.

For the material constitutive model, we use isotropic Neo-Hookean material as
detailed in [33]. The mechanical properties $\Theta$, including Young’s
modulus $E$, Poisson’s ratio $\nu$, and mass density $\rho$, are set by users.
These values can be specified directly through numerical input or chosen from
a collection of pre-established material options, such as plastic or rubber.
We consider gravity and fixed attachment forces as options for external
forces. Gravity is always included to reflect its omnipresence in the real
world. The use of fixed attachment forces depends on the specific needs of the
application, for instance, anchoring an object at a designated site. Detailed
formulations for both force types are provided in Appendix F.

##  3 Evaluation

In this section, we present evidence that our approach enhances the physical
compatibility of 3D objects produced using state-of-the-art single-view
reconstruction techniques. We conduct a series of quantitative evaluations
using five metrics (Sec. 3.1) to compare the physical compatibility of shapes
optimized by our framework against those produced by existing methods without
our method (Sec. 3.2). We also provide qualitative comparisons to demonstrate
to the effectiveness of our approach (Sec. 3.3). Furthermore, we explore the
practical applications of our method by illustrating how it enables the
reconstruction of diverse 3D shapes with different material properties from
the same single image, and by demonstrating that our optimized shapes are
readily adaptable for dynamic simulations and fabrication (Sec. 3.4).

###  3.1 Baselines and Evaluation Protocol

Existing metrics for evaluating single-view reconstruction methods primarily
focus on the visual appearance of the objects. Measures such as PSNR and SSIM
are used to assess image fidelity, while chamfer distance and volume IoU
evaluate geometric quality. However, these metrics do not consider the
underlying physics principles that govern the behavior of 3D objects.
Consequently, they are insufficient for evaluating the physical compatibility
of reconstructed shapes, a crucial aspect for applications requiring accurate
physical interactions and structural stability.

Table 1:  Quantitative results on four metrics evaluating physical
compatibility. We apply our pipeline to five single-image reconstruction
techniques and assess our metrics on both the initial shapes from these
methods (Baseline) and the optimized shapes from the integration of our
framework with each baseline (Ours). Our method demonstrates quantitative
improvements in mean stress, stability rate, and image fidelity across all
benchmarks. Among all methods, TetSphere integrated with our framework
achieves superior performance across all evaluation metrics. This can be
attributed to the explicit volumetric representation used in TetSphere. The
mean and standard deviation are calculated across all examples for each
method. A higher deviation in Mean Stress suggests a larger variance in
structural thickness and curvature, while a higher deviation in Img. Loss
indicates a larger variance in static shape deformation.

Method | Init. Geo. | $\\#$CC. $\downarrow$ |  Mean Stress $\downarrow$ (kPa) |  Standable. $\uparrow$ (%) |  Img. Loss $\downarrow$
---|---|---|---|---|---
Wonder3D | Baseline | NeuS | 2.54 $\pm$ 2.64 |  10.68 $\pm$ 17.47 | 6.9 |  0.073 $\pm$ 0.063
Ours |  0.45 $\pm$ 0.96 | 72.4 |  0.069 $\pm$ 0.048
LGM | Baseline |  Gaussian splatting | 2.67 $\pm$ 2.13 |  1.14 $\pm$ 2.03 | 20.3 |  0.121 $\pm$ 0.091
Ours |  1.01 $\pm$ 1.34 | 85.9 |  0.116 $\pm$ 0.065
MeshLRM | Baseline |  surface mesh | 1.55$\pm$ 2.13 |  0.54 $\pm$ 1.41 | 28.6 |  0.065 $\pm$ 0.042
Ours |  0.38 $\pm$ 1.05 | 73.8 |  0.064 $\pm$ 0.042
TripoSR | Baseline | NeRF |  1.43 $\pm$ 1.12 |  0.29 $\pm$ 1.28 | 24.2 |  0.066 $\pm$ 0.047
Ours |  0.22 $\pm$ 0.94 | 80.6 |  0.059 $\pm$ 0.039
TetSphere | Baseline | tet-sphere | 1.00 $\pm$ 0.00 |  0.22 $\pm$ 0.51 | 32.8 |  0.061 $\pm$ 0.045
Ours | 0.19 $\pm$ 0.78 | 92.4 | 0.057 $\pm$ 0.040

Figure 3:  Quantitative results on fracture rate. We plot the relationship
between the fracture rate and the maximum stress threshold across five single-
image reconstruction methods. The shapes optimized with our framework exhibit
a consistently lower fracture rate compared to those shapes obtained without
our pipeline. MeshLRM and TripoSR feature prevalent thin structures in their
reconstructed shapes, whereas our approach significantly reduces the fracture
rate in both cases.

#### Metrics.

To address this oversight, we draw inspiration from the field of finite
element analysis [2] and introduce five novel metrics specifically designed to
assess the physical compatibility of 3D models comprehensively. These metrics
are tailored to ensure a more thorough evaluation of method performance in
real-world scenarios with rich physics:

  * •

Number of Connected Components ($\\#$CC.) evaluates the structural integrity
of the object. Physical objects should not have floating or disconnected
structures, ideally consisting of one single connected component.

  * •

Mean Stress calculates the average von Mises stress [26] across all tetrahedra
of all objects. It measures the extent of physical deformation. Under the same
external interactions, higher mean stress indicates a greater likelihood of
fracture and the existence of unrealistic thin structures.

  * •

Percentage of Standability (Standable.) assesses whether the object can
maintain stability under gravity, remaining upright without toppling. A
standable object is one that effectively supports itself against gravitational
forces.

  * •

Matching loss (Img. Loss) calculates the $l_{1}$ difference between the
rendered image of the object after applying gravity and the input target
image, quantifying the deviation of the physical object from the desired shape
due to physical influences.

  * •

Fracture Rate measures the number of tetrahedral elements that exceed a
predefined stress threshold, potentially leading to fractures. The resilience
of a method against physical stresses is quantified using a degradation curve,
with more physically reliable methods exhibiting a smaller area under the
curve for the fracture rate.

Figure 4:  Qualitative results on physical compatibility optimization. Left:
Rest shapes optimized using our approach result in static shapes that closely
match the input images when subjected to gravity. In contrast, shapes without
the optimization fail to replicate the geometry in the input image. Right: our
optimization process ensures that the optimized shapes are capable of
supporting themselves, whereas the baseline methods fail to achieve this
stability.

Baselines. We consider five single-view reconstruction baselines in our
evaluation, each associated with a distinct geometry representation: Wonder3D
[23] with NeuS, LGM [38] with Gaussian splatting, MeshLRM [42] with surface
mesh, TripoSR [39] with NeRF triplane, and TetSphere [11] with tetrahedral
spheres. For the baseline results, we used the publicly available inference
code to reconstruct the 3D objects.222For MeshLRM, since the pre-trained model
is not publicly available yet, we obtained the reconstructed shapes directly
from the authors for use in our study. To demonstrate the versatility of our
method, we integrated our physical compatibility optimization framework with
all five baseline models and reported the results to ensure a fair comparison.
The implementation details of our framework are provided in Appendix D.

Evaluation Datasets. The evaluation dataset was sourced from Objaverse [9]. We
initially randomly selected approximately $200$ shapes from the categories of
plants, animals, and characters – categories that demand greater physical
compatibility. Single-view images were rendered using the publicly released
code by the authors of Objaverse333https://github.com/allenai/objaverse-
rendering. Subsequently, these images were used to reconstruct 3D objects
using the baseline methods mentioned earlier. We filtered out shapes of
extremely poor quality, specifically those with more than $8$ connected
components. This process resulted in a final set of $100$ shapes for detailed
evaluation.

Despite these shapes being a part of the training data for most baseline
methods, our evaluation focuses on assessing the physical compatibility – a
factor overlooked by these methods. The results obtained from this dataset
provide valuable insights and observations on the physical compatibility of
each method, demonstrating the practical effectiveness of our approach.

Figure 5:  Ablation study on Young’s modulus. By changing the material
properties, our method can produce various rest-shape geometries (top), which
all result in the same static shapes that match the input image (middle).
Although these static shapes appear identical under static equilibrium, they
exhibit different deformation when subjected to the same compression forces
exerted by the yellow block, attributable to the differences in their material
properties (bottom).

###  3.2 Quantitative Results

Table 1 shows the quantitative results for four out of five metrics evaluated
for both baselines and those integrated with our physical compatibility
optimization. Fig. 3 shows the curve of fracture rate.

Our quantitative analysis yields several observations: 1) The underlying
geometry representation significantly impacts the structural integrity of
reconstructed shapes, as evidenced by the number of connected components
($\\#$CC.). LGM, using a point cloud representation, exhibits the poorest
structural integrity, often resulting in floating structures due to its
inability to differentiate the interior from the exterior of a 3D object. In
contrast, TetSphere, with its volumetric representation, maintains the most
integral structure. 2) Both MeshLRM and TripoSR generally produce more
physically stable 3D objects, as indicated by Mean Stress and Standability
(Standable.) metrics. However, they tend to diverge under gravity, as shown by
the Matching Loss metric (Img. Loss), compared to TetSphere. 3) Notably, our
method consistently enhances the physical compatibility performance across all
baselines. The improvement is particularly significant for Wonder3D and
MeshLRM. Wonder3D typically generates multi-view images before reconstructing
the 3D shape, which can lead to thin structures due to inconsistencies across
the views. Similarly, MeshLRM’s reliance on surface mesh could often result in
thin structures. Our method strengthens the physical robustness for both
cases. 4) Our method also enhances the structure robustness to fracture, as
demonstrated in Fig. 3. It notably improves the performance of both MeshLRM
and TripoSR in reducing fracture rates.

Figure 6:  Applications of physically compatible objects. Left: Our optimized
physical objects is simulation-ready and can be seamlessly integrated into
dynamic simulation pipeline to produce complex dynamics and motions. Right:
Real-world validation using 3D printing shows that shapes optimized using our
method closely replicate the input images, demonstrating the practical
effectiveness of our method in manufacturing.

###  3.3 Qualitative Results

Fig. 4 and more qualitative results in Appendix 4 illustrate the effectiveness
of our physical compatibility optimization. Without optimization, the static
shapes behave undesirably under general physical principles: they either sag
excessively under gravity, diverging from the geometry depicted in the input
image, or fail to remain upright, toppling over. Our optimization method
incorporates physical principles to ensure that the optimized rest shapes are
self-supporting and stable, and match the input images under static
equilibrium.

###  3.4 Analysis

#### Ablation study on Young’s Modulus.

We investigate the influence of predefined mechanical material properties,
particularly Young’s modulus, on the optimized rest shapes and their physical
behaviors. Using the same input image, we obtained six optimized rest shapes
with varying Young’s modulus values within our framework with TetSphere. As
shown in Fig. 5, although the optimized rest-shape geometries vary, they all
deform to the same static geometry under the influence of gravity, matching
the input image. Moreover, the physical responses to identical external
forces, such as compression by a box, differ due to the variations in material
properties. These results highlight how the explicit decomposition of physical
attributes in our framework expands the controllability of object modeling,
allowing for diverse physical behaviors under uniform external forces.

#### Application to dynamic simulation.

The immediate output of our method is a simulation-ready rest-shape geometry,
which can be seamlessly integrated into a simulation pipeline to produce
complex dynamics and motions. Fig. 6 (left) and the accompanying video in the
Supplementary Material illustrate three plants modeled using our framework,
demonstrating their behavior under gravity and complex interactions.
Implementation details of this simulation are provided in Appendix F. These
examples underscore the practical utility of our method for generating
physically realistic dynamics and simulations.

#### Application to fabrication.

We further evaluate our method in real world by fabricating three shapes using
3D printing, both with and without optimization. The results, shown in Fig. 6
(right), with detailed implementation procedures available in Appendix E,
demonstrate that the 3D printed shapes align with our computational results.
These real-world experiments demonstrate the practical effectiveness and
validate the physical realism of the objects produced by our method.

##  4 Conclusion

In this work, we introduced physical compatibility optimization for
reconstructing a physical object from a single image. Our method decomposes
three orthogonal attributes governing physical behavior: mechanical
properties, external forces, and rest-shape geometry. Unlike existing methods
that often ignore one or more dimensions, our framework holistically considers
all three factors, allowing for diverse rest-shape geometries from the same
input image by varying object stiffness and external forces. We formulate
physical compatibility optimization as a constrained optimization problem by
integrating static equilibrium as a hard constraint. Our approach produces
physical objects that match the geometry depicted in the input image under
external forces and remain stable under gravity. Both quantitative and
qualitative evaluations demonstrated improvements in physical compatibility
over existing baselines. Our method’s versatility is evident through its
integration with various single-view reconstruction methods and its practical
applications in dynamic simulations and 3D printing.

## References

  * Aifantis [1987] E. C. Aifantis.  The physics of plastic deformation.  _International journal of plasticity_ , 3(3):211–247, 1987.
  * Allaire [2007] G. Allaire.  _Numerical analysis and optimization: an introduction to mathematical modelling and numerical simulation_.  OUP Oxford, 2007.
  * Bell et al. [2014] S. Bell, K. Bala, and N. Snavely.  Intrinsic images in the wild.  _ACM Transactions on Graphics (TOG)_ , 33(4):1–12, 2014.
  * Bermano et al. [2017] A. H. Bermano, T. Funkhouser, and S. Rusinkiewicz.  State of the art in methods and representations for fabrication-aware design.  In _Computer Graphics Forum_ , volume 36, pages 509–535. Wiley Online Library, 2017.
  * Botsch and Sorkine [2007] M. Botsch and O. Sorkine.  On linear variational surface deformation methods.  _IEEE transactions on visualization and computer graphics_ , 14(1):213–230, 2007.
  * Burczyński et al. [1997] T. Burczyński, J. Kane, and C. Balakrishna.  Comparison of shape design sensitivity analysis formulations via material derivative-adjoint variable and implicit differentiation techniques for 3-d and 2-d curved boundary element.  _Computer methods in applied mechanics and engineering_ , 142(1-2):89–109, 1997.
  * Charatan et al. [2023] D. Charatan, S. Li, A. Tagliasacchi, and V. Sitzmann.  pixelsplat: 3d gaussian splats from image pairs for scalable generalizable 3d reconstruction.  _arXiv preprint arXiv:2312.12337_ , 2023.
  * Chen et al. [2014] X. Chen, C. Zheng, W. Xu, and K. Zhou.  An asymptotic numerical method for inverse elastic shape design.  _ACM Transactions on Graphics (TOG)_ , 33(4):1–11, 2014.
  * Deitke et al. [2023] M. Deitke, D. Schwenk, J. Salvador, L. Weihs, O. Michel, E. VanderBilt, L. Schmidt, K. Ehsani, A. Kembhavi, and A. Farhadi.  Objaverse: A universe of annotated 3d objects.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 13142–13153, 2023.
  * Feng et al. [2023] Y. Feng, Y. Shang, X. Li, T. Shao, C. Jiang, and Y. Yang.  Pie-nerf: Physics-based interactive elastodynamics with nerf.  _arXiv preprint arXiv:2311.13099_ , 2023.
  * [11] M. Guo, B. Wang, K. He, and W. Matusik.  Tetsphere splatting: Representing high-quality geometry with lagrangian volumetric meshes.
  * He and Wang [2023] Z. He and T. Wang.  Openlrm: Open-source large reconstruction models.  https://github.com/3DTopia/OpenLRM, 2023.
  * Hong et al. [2023] Y. Hong, K. Zhang, J. Gu, S. Bi, Y. Zhou, D. Liu, F. Liu, K. Sunkavalli, T. Bui, and H. Tan.  LRM: Large reconstruction model for single image to 3D.  Nov. 2023.
  * Hsu et al. [2022] J. Hsu, N. Truong, C. Yuksel, and K. Wu.  A general two-stage initialization for sag-free deformable simulations.  _ACM Transactions on Graphics (TOG)_ , 41(4):1–13, 2022.
  * Hu et al. [2018] Y. Hu, Q. Zhou, X. Gao, A. Jacobson, D. Zorin, and D. Panozzo.  Tetrahedral meshing in the wild.  _ACM Trans. Graph._ , 37(4):60:1–60:14, July 2018.  ISSN 0730-0301.  doi: 10.1145/3197517.3201353.  URL http://doi.acm.org/10.1145/3197517.3201353.
  * Kerbl et al. [2023] B. Kerbl, G. Kopanas, T. Leimkühler, and G. Drettakis.  3d gaussian splatting for real-time radiance field rendering.  _ACM Transactions on Graphics_ , 42(4), 2023.
  * Kingma and Ba [2015] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _ICLR_ , 2015.
  * Li et al. [2020] M. Li, Z. Ferguson, T. Schneider, T. R. Langlois, D. Zorin, D. Panozzo, C. Jiang, and D. M. Kaufman.  Incremental potential contact: intersection-and inversion-free, large-deformation dynamics.  _ACM Trans. Graph._ , 39(4):49, 2020.
  * Li et al. [2022] X. Li, Y.-L. Qiao, P. Y. Chen, K. M. Jatavallabhula, M. Lin, C. Jiang, and C. Gan.  Pac-nerf: Physics augmented continuum neural radiance fields for geometry-agnostic system identification.  2022\.
  * Liu et al. [2023a] M. Liu, C. Xu, H. Jin, L. Chen, V. T. Mukund, Z. Xu, and H. Su.  One-2-3-45: Any single image to 3D mesh in 45 seconds without Per-Shape optimization.  June 2023a.
  * Liu et al. [2023b] R. Liu, R. Wu, B. V. Hoorick, P. Tokmakov, S. Zakharov, and C. Vondrick.  Zero-1-to-3: Zero-shot one image to 3d object, 2023b.
  * Liu et al. [2023c] Y. Liu, C. Lin, Z. Zeng, X. Long, L. Liu, T. Komura, and W. Wang.  SyncDreamer: Generating multiview-consistent images from a single-view image.  Sept. 2023c.
  * Long et al. [2023] X. Long, Y.-C. Guo, C. Lin, Y. Liu, Z. Dou, L. Liu, Y. Ma, S.-H. Zhang, M. Habermann, C. Theobalt, and W. Wang.  Wonder3D: Single image to 3D using Cross-Domain diffusion.  Oct. 2023.
  * Melas-Kyriazi et al. [2023] L. Melas-Kyriazi, C. Rupprecht, I. Laina, and A. Vedaldi.  RealFusion: 360° reconstruction of any object from a single image.  Feb. 2023.
  * Mildenhall et al. [2020] B. Mildenhall, P. P. Srinivasan, M. Tancik, J. T. Barron, R. Ramamoorthi, and R. Ng.  Nerf: Representing scenes as neural radiance fields for view synthesis.  _ECCV_ , 2020.
  * Mises [1913] R. v. Mises.  Mechanik der festen körper im plastisch-deformablen zustand.  _Nachrichten von der Gesellschaft der Wissenschaften zu Göttingen, Mathematisch-Physikalische Klasse_ , 1913:582–592, 1913.
  * Nicolet et al. [2021] B. Nicolet, A. Jacobson, and W. Jakob.  Large steps in inverse rendering of geometry.  _ACM Transactions on Graphics (TOG)_ , 40(6):1–13, 2021.
  * Perlin [2002] K. Perlin.  Improving noise.  In _Proceedings of the 29th annual conference on Computer graphics and interactive techniques_ , pages 681–682, 2002.
  * Poole et al. [2022] B. Poole, A. Jain, J. T. Barron, and B. Mildenhall.  Dreamfusion: Text-to-3d using 2d diffusion.  _arXiv preprint arXiv:2209.14988_ , 2022.
  * Qian et al. [2023] G. Qian, J. Mai, A. Hamdi, J. Ren, A. Siarohin, B. Li, H.-Y. Lee, I. Skorokhodov, P. Wonka, S. Tulyakov, and B. Ghanem.  Magic123: One image to high-quality 3d object generation using both 2d and 3d diffusion priors.  _arXiv preprint arXiv:2306.17843_ , 2023.
  * Shue et al. [2023] J. R. Shue, E. R. Chan, R. Po, Z. Ankner, J. Wu, and G. Wetzstein.  3d neural field generation using triplane diffusion.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 20875–20886, 2023.
  * Sifakis and Barbic [2012] E. Sifakis and J. Barbic.  Fem simulation of 3d deformable solids: a practitioner’s guide to theory, discretization and model reduction.  In _Acm siggraph 2012 courses_ , pages 1–50. 2012.
  * Smith et al. [2018] B. Smith, F. D. Goes, and T. Kim.  Stable neo-hookean flesh simulation.  _ACM Transactions on Graphics (TOG)_ , 37(2):1–15, 2018.
  * Smith and Schaefer [2015] J. Smith and S. Schaefer.  Bijective parameterization with free boundaries.  _ACM Transactions on Graphics (TOG)_ , 34(4):1–9, 2015.
  * Standley et al. [2017] T. Standley, O. Sener, D. Chen, and S. Savarese.  image2mass: Estimating the mass of an object from its image.  In _Conference on Robot Learning_ , pages 324–333. PMLR, 2017.
  * Szymanowicz et al. [2024] S. Szymanowicz, C. Rupprecht, and A. Vedaldi.  Splatter image: Ultra-fast single-view 3d reconstruction.  2024\.
  * Tang et al. [2023] J. Tang, J. Ren, H. Zhou, Z. Liu, and G. Zeng.  DreamGaussian: Generative gaussian splatting for efficient 3D content creation.  Sept. 2023.
  * Tang et al. [2024] J. Tang, Z. Chen, X. Chen, T. Wang, G. Zeng, and Z. Liu.  LGM: Large Multi-View gaussian model for High-Resolution 3D content creation.  Feb. 2024.
  * Tochilkin et al. [2024] D. Tochilkin, D. Pankratz, Z. Liu, Z. Huang, A. Letts, Y. Li, D. Liang, C. Laforte, V. Jampani, and Y.-P. Cao.  Triposr: Fast 3d object reconstruction from a single image.  _arXiv preprint arXiv:2403.02151_ , 2024.
  * Wang et al. [2021a] B. Wang, G. Matcuk, and J. Barbič.  Modeling of personalized anatomy using plastic strains.  _ACM Transactions on Graphics (TOG)_ , 40(2):1–21, 2021a.
  * Wang et al. [2021b] P. Wang, L. Liu, Y. Liu, C. Theobalt, T. Komura, and W. Wang.  NeuS: Learning neural implicit surfaces by volume rendering for multi-view reconstruction.  June 2021b.
  * Wei et al. [2024] X. Wei, K. Zhang, S. Bi, H. Tan, F. Luan, V. Deschaintre, K. Sunkavalli, H. Su, and Z. Xu.  Meshlrm: Large reconstruction model for high-quality mesh.  _arXiv preprint arXiv:2404.12385_ , 2024.
  * Wu et al. [2023] R. Wu, B. Mildenhall, P. Henzler, K. Park, R. Gao, D. Watson, P. P. Srinivasan, D. Verbin, J. T. Barron, B. Poole, et al.  Reconfusion: 3d reconstruction with diffusion priors.  _arXiv preprint arXiv:2312.02981_ , 2023.
  * Xie et al. [2023] T. Xie, Z. Zong, Y. Qiu, X. Li, Y. Feng, Y. Yang, and C. Jiang.  Physgaussian: Physics-integrated 3d gaussians for generative dynamics.  _arXiv preprint arXiv:2311.12198_ , 2023.
  * Xu et al. [2024] Y. Xu, Z. Shi, W. Yifan, H. Chen, C. Yang, S. Peng, Y. Shen, and G. Wetzstein.  Grm: Large gaussian reconstruction model for efficient 3d reconstruction and generation.  _arXiv preprint arXiv:2403.14621_ , 2024.
  * Yu et al. [2021] A. Yu, V. Ye, M. Tancik, and A. Kanazawa.  pixelNeRF: Neural radiance fields from one or few images.  2021\.
  * Zhai et al. [2024] A. J. Zhai, Y. Shen, E. Y. Chen, G. X. Wang, X. Wang, S. Wang, K. Guan, and S. Wang.  Physical property understanding from language-embedded feature fields.  _arXiv preprint arXiv:2404.04242_ , 2024.
  * Zhang et al. [2024] K. Zhang, S. Bi, H. Tan, Y. Xiangli, N. Zhao, K. Sunkavalli, and Z. Xu.  Gs-lrm: Large reconstruction model for 3d gaussian splatting.  _arXiv preprint arXiv:2404.19702_ , 2024.
  * Zhong et al. [2024] L. Zhong, H.-X. Yu, J. Wu, and Y. Li.  Reconstruction and simulation of elastic objects with spring-mass 3d gaussians.  _arXiv preprint arXiv:2403.09434_ , 2024.

##  Appendix A Additional Qualitative Results

Figure 7: Additional qualitative results of physical compatibility
optimization (part 1/2).

Figure 8: Additional qualitative results of physical compatibility
optimization (part 2/2).

Figure 7 and 8 show additional results of our physical compatibility
optimization.

##  Appendix B Plastic Strain Field $\mathbf{F_{p}}$

To enhance the understanding of our framework without compromising
generalizability, let us consider $\mathcal{M}_{\mathrm{init}}$ to be a
tetrahedral mesh composed of a single element and four vertices. When subject
to static equilibrium influenced by gravity, the object adheres to the
equation:

$\mathbf{f}_{\mathrm{int}}(\mathbf{x},\phi(\mathbf{F_{p}};\mathbf{X}_{\mathrm{init}});\Theta)=\mathbf{M}\mathbf{g},$
(6)

where $\mathbf{f}_{\mathrm{int}}(\cdot,\cdot)$ denotes the elastic force
(internal force), $\mathbf{M}$ is the mass matrix, and $\mathbf{g}$ denotes
the gravity acceleration. To compute this force, we first consider the elastic
energy $\mathcal{E}$. The definition of elastic energy unfolds as follows:

$\displaystyle\mathcal{E}(\mathbf{F_{e}},\mathbf{F_{p}};\Theta)$
$\displaystyle=V(\mathbf{F_{p}})\Phi(\mathbf{F_{e}};\Theta),$ $\displaystyle
V(\mathbf{F_{p}})$
$\displaystyle=V_{\mathrm{init}}\mathrm{det}(\mathbf{F_{p}}),$
$\displaystyle\mathbf{F_{e}}$ $\displaystyle=\mathbf{F}\mathbf{F_{p}}^{-1},$
$\displaystyle\mathbf{F}$
$\displaystyle=\partial\mathbf{x}/\partial\mathbf{X}_{\mathrm{init}},$

where $V(\mathbf{F_{p}})$ represents the volume of the element under plastic
strain, $V_{\mathrm{init}}$ is the initial volume of the element,
$\mathbf{F_{e}}$ denotes the elastic deformation gradient, $\mathbf{F}$
represent the total deformation gradient, and $\Phi(\cdot;\Theta)$ the elastic
energy density function. This deformation gradient $\mathbf{F}$ is computed
through standard methodology [32].

Consequently, the derivation of the elastic force encapsulates the computation
of the first-order partial derivative of the elastic energy with respect to
the vertex positions:

$\displaystyle\mathbf{f}_{\mathrm{int}}(\mathbf{x},\phi(\mathbf{F_{p}};\mathbf{X}_{\mathrm{init}});\Theta)$
$\displaystyle:=\frac{\partial\mathcal{E}(\mathbf{F_{e}}(\mathbf{x}),\mathbf{F_{p}};\Theta)}{\partial\mathbf{x}}$
$\displaystyle=V(\mathbf{F_{p}})\frac{\partial\Phi}{\partial\mathbf{F_{e}}}\colon\frac{\partial\mathbf{F}}{\partial\mathbf{x}}\mathbf{F_{p}}^{-1}.$

Notably, given the linear dependence of $\mathbf{F}$ on $\mathbf{x}$,
$\frac{\partial\mathbf{F}}{\partial\mathbf{x}}$ remains constant.

Given $\mathbf{F_{p}}$ and $\mathbf{X}_{\mathrm{init}}$ as inputs, the
solution to Eq. 6 is the static shape,
$\mathbf{x}=\mathbf{x}_{\mathrm{static}}$. Likewise, to calculate
$\mathbf{X}_{\mathrm{rest}}$ from $\mathbf{F_{p}}$ and
$\mathbf{X}_{\mathrm{init}}$ in Eq. 3, we solve a similar equation with zero
external force.

$\displaystyle\mathbf{f}_{\mathrm{int}}(\mathbf{x},\phi(\mathbf{F_{p}};\mathbf{X}_{\mathrm{init}});\Theta)$
$\displaystyle=\mathbf{0},$

where the solution to this equation is $\mathbf{X}_{\mathrm{rest}}$.

Considering the elastic energy, the translation of
$\mathbf{X}_{\mathrm{init}}$ does not alter the deformation gradient
$\mathbf{F}$. Consequently, $\mathbf{F_{p}}$ remain unaffected and exhibit
translation invariance. In terms of the elastic force, it maintains
translation invariance as well, since $\mathbf{F}$ is not affected by any
shift in $\mathbf{X}_{\mathrm{init}}$.

Finally, by using isotropic materials, our approach enables a further
reduction in the DOFs of $\mathbf{F_{p}}$. Let us denote $\mathbf{F_{p}}$ as
$\mathbf{F_{p}}=\mathbf{R}\mathbf{S}$. The elastic deformation gradient is
then derived as
$\mathbf{F_{e}}=\mathbf{F}(\mathbf{R}\mathbf{S})^{-1}=\mathbf{F}\mathbf{S}^{-1}\mathbf{R}^{-1}$.
Given the invariance property
$\Phi(\mathbf{F_{e}};\theta)=\Phi(\mathbf{F_{e}}\mathbf{R};\theta)$, which
constantly holds for isotropic materials, the rotation component $\mathbf{R}$
becomes redundant and can be excluded from the formulation. This
simplification implies that the only requirement for $\mathbf{F_{p}}$ is to be
a symmetric matrix. During the optimization process, this property facilitates
the prevention of the inversion: In order to ensure that
$\mathrm{det}(\mathbf{F_{p}})>0$, we can simply adjust the eigenvalues of
$\mathbf{F_{p}}$ to make they remain positive. This adjustment is crucial for
the rest mesh $\mathbf{X}_{\mathrm{rest}}$ to maintain in the non-inverted
state.

##  Appendix C Computation of Gradient

By differentiating the constraint in Eq. 2.2 with respect to $\mathbf{F_{p}}$,
we obtain

$\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{F_{p}}}+\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{x}_{\mathrm{static}}}\frac{\partial\mathbf{x}_{\mathrm{static}}}{\partial\mathbf{F_{p}}}=0.$
(7)

Then, we have

$\frac{\partial\mathbf{x}_{\mathrm{static}}}{\partial\mathbf{F_{p}}}=-[\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{x}_{\mathrm{static}}}]^{-1}\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{F_{p}}}.$
(8)

Substituting the result into the objective in Eq. 2.2, we get

$\displaystyle\frac{\partial\mathcal{J}}{\partial\mathbf{F_{p}}}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\mathbf{F_{p}}}+\frac{\partial\mathcal{L}_{\mathrm{reg}}}{\partial\mathbf{F_{p}}}$
$\displaystyle=\frac{\partial\mathcal{L}}{\partial\mathbf{x}_{\mathrm{static}}}\frac{\partial\mathbf{x}_{\mathrm{static}}}{\partial\mathbf{F_{p}}}+\frac{\partial\mathcal{L}_{\mathrm{reg}}}{\partial\mathbf{F_{p}}}$
$\displaystyle=-\frac{\partial\mathcal{L}}{\partial\mathbf{x}_{\mathrm{static}}}[\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{x}_{\mathrm{static}}}]^{-1}\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{F_{p}}}+\frac{\partial\mathcal{L}_{\mathrm{reg}}}{\partial\mathbf{F_{p}}},$
(9)

which is the gradient with respect to $\mathbf{F_{p}}$. In practice,
$\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{x}_{\mathrm{static}}}$
and $\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{F_{p}}}$ are
stored as sparse matrices and computed based on [40]. Considering about the
performance, we first compute
$\frac{\partial\mathcal{L}}{\partial\mathbf{x}_{\mathrm{static}}}[\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{x}_{\mathrm{static}}}]^{-1}$
using sparse linear solver. This results in a dense vector with size $3N$. We
then multiply it with
$\frac{\partial\mathbf{f}_{\mathrm{net}}}{\partial\mathbf{F_{p}}}$.

##  Appendix D Implementation Details of Evaluation

To evaluate the physical compatibility of baseline methods, which often
produce shapes comprising multiple connected components, we first extract the
largest connected component from each mesh. All meshes are then normalized to
the unit cube. Notably, the reconstructed shapes from TripoSR and Wonder3D are
not axis-aligned; thus, we manually rotate these shapes to ensure the head
points towards the $z$-axis in the world coordinate space. For integrating our
physical compatibility framework, We use two sets of Young’s modulus,
$E=5\times 10^{4}\mathrm{Pa}$ and $E=5\times 10^{5}\mathrm{Pa}$, which are
selected based on whether the shape would become overly soft, potentially
leading to static equilibrium failure due to excessive stress causing
numerical bounds to be exceeded. Poisson’s ratio $\nu=0.45$ and mass density
$\rho=1000\mathrm{kg/m^{3}}$ are consistent across all meshes. Evaluation
metrics require solving for static equilibrium Eq. 1. We employ the Newton-
Raphson solver with line search, setting the maximum number of iterations to
be $200$. For optimizing Eq. 2.2, we use gradient descent and allow up to
$1000$ iterations. Our experiments run on a desktop PC with an AMD Ryzen 9
5950X 16-core CPU and 64GB RAM. The average runtime for this optimization
process is approximately $80$ seconds.

##  Appendix E Implementation Details of 3D Printing

The selected model shapes were 3D printed using stereolithography (Form3;
Formlabs, $100$ $\mu$m layer thickness) to create the flexible designs (using
Flexible 80A, tensile modulus ${<}3$ MPa, 100% strain to failure) and rigid
designs (using White Resin V4; tensile modulus $1.6$ GPa), both without post-
curing. The flexible flowers are $55$ and $65$ mm in height and the rigid
goose is $50$ mm in length. Shapes with and without optimization were printed
with similar support structures designed to preserve delicate features.

##  Appendix F Dynamic Simulation of Deformable Objects

We model each solid deformable object using FEM with hyperelastic materials
for dynamic simulation. Then, we solve the standard partial differential
equation (PDE) for dynamic FEM simulation:

$M\ddot{x}+D(x)\dot{x}+f_{\mathrm{elastic}}(x)+f_{\mathrm{attachment}}(x)+f_{\mathrm{contact}}(x)=Mg,$
(10)

where $x$ represents the node positions within the finite element meshes – we
use tetrahedral meshes – of the objects, $M$ denotes the mass matrix, $D$ is
the Rayleigh damping matrix, $f_{\mathrm{elastic}}(\cdot)$ is the hyperelastic
forces, $f_{\mathrm{attachment}}(x)$ is the attachment forces that constrain
the objects to a specific location, and $f_{\mathrm{contact}}(\cdot)$ denotes
the contact forces between surfaces. We employ the implicit backward Euler
method for time discretization, transforming the PDE into:

$A^{n}x^{n+1}+b^{n}+f_{\mathrm{elastic}}(x^{n+1})+f_{\mathrm{attachment}}(x^{n+1})+f_{\mathrm{contact}}(x^{n+1})=0,$
(11)

where $x^{n+1}$ is the position vector at timestep $(n+1)$, $A^{n}$ and
$b^{n}$ is a constant matrix and vector, respectively, derived from values at
timestep $n$, Finally, we solve this nonlinear equation using Newton’s method
at each timestep.

The hyperelastic material selected for the deformable objects is the same as
the one used for the rest shape optimization [33] in Sec. 2. Attachment forces
are modeled as spring forces
$f_{\mathrm{attachment}}(x)=k_{a}(Sx-\bar{x}(t))$, where $k_{a}$ is the
stiffness of the spring, the selection matrix $S$ selects the attached
vertices, and $\bar{x}(t)$ denotes the target attachment locations at time
$t$. Contact forces are generated from penalizing any vertex penetration into
the contact surface, expressed as $f=k_{c}d$, where $k_{c}$ represents the
contact stiffness and $d$ denotes the penetration depth, with $d=0$ in the
absence of contact. This gives the normal contact forces. Friction forces are
computed following the methods outlined in [18]. Then, the total contact force
$f_{\mathrm{contact}}$ is the sum of normal contact forces and friction
forces.

For the dynamic simulation in Figure 7, the attachment of each plant is
defined as the bottom part of each pot. We keyframe-animate the trajectory of
attachment vertices $\bar{x}(t)$. Gravity is enabled throughout the entire
simulation. At the end of the sequence, we apply wind forces to the plants,
computed using 4D Perlin Noise [28].

##  Appendix G Related Work

#### Single-view 3D reconstruction.

Recent strides in single-view 3D reconstruction have mainly been fueled by
data-driven methods, paralleled by advancements in 3D geometry representation,
including NeRF [25], NeuS [41], triplanes [31], Gaussian splatting [16],
surface meshes [27], and tet-spheres [11]. These developments have
significantly enhanced the geometric quality and visual fidelity of the
reconstructed 3D shapes. There are primarily two types of single-view
reconstruction methods: 1) Test-time optimization-based methods [29, 23, 37,
43], use multiview diffusion models [21] and iteratively reconstruct 3D scenes
using these diffusion priors. 2) Feedforward methods [13, 46, 36, 7, 42, 48]
leverage large datasets and learn general 3D priors for shape reconstruction
to enable efficient one-step 3D reconstruction from single or sparse views.
Unlike the aforementioned methods, our work emphasizes the integration of
physical modeling into the reconstruction process. This integration
distinguishes our work by ensuring that the resulting 3D shapes are not only
visually accurate but also physically plausible under real-world conditions.

#### Physics-based 3D modeling.

There has been an increasing interest in incorporating physics into 3D shape
modeling. While many approaches utilize video input, which offers a richer
temporal context for inferring physical properties such as material parameters
[49] and geometry [19], others approach the problem by first reconstructing an
object’s geometry from multi-view images and subsequently applying physical
simulations [10, 44]. Additionally, several studies have explored extracting
physical information from static images [47, 3, 35], using data-driven
techniques to estimate properties like shading, mass, and material. In
contrast, our work incorporates physical principles, specifically static
equilibrium, as hard constraints within the reconstruction process. This
integration allows for the optimization of 3D models that adhere to desired
physical behaviors depicted by the image.

#### Fabrication-aware shape design.

Originating from the computer graphics community, fabrication-aware shape
design systems enable designers to specify higher-level objectives – such as
structural integrity, deformation, and appearance – with the final shape as
the output of the computational system [4]. Related methodologies in this
domain, particularly those addressing static equilibrium, include inverse
elastic shape design [8] and sag-free initialization [14]. However, these
approaches typically require a manually created initial geometry, whereas our
work aims to construct the physical object directly from a single input image.

##  Appendix H Limitations and Future Work

One limitation of our framework is its reliance on predefined material
properties and external forces as inputs. Although this provides
controllability of the final optimized rest-shape geometry, automating the
extraction of these parameters from a single image presents a potential avenue
for future work. Moreover, our method relies on the use of a tetrahedral mesh,
which is derived by tetrahedralizing the output geometry produced by baseline
methods. A natural extension of our work is the development of a
differentiable converter that can transform any geometric representation into
a tetrahedral mesh. This would enable future research where our physical
compatibility optimization could be integrated into a pre-trained large
reconstruction model, which could then be fine-tuned to directly produce
physically compatible 3D objects. Lastly, our current methodology focuses
solely on physical objects in a state of static equilibrium. Exploring the
reconstruction of 3D objects undergoing dynamics captured from video is an
intriguing prospect for future research.

##  Appendix I Broader Impacts

Our research presents a computational framework for reconstructing physical
objects from single images. This advancement holds significant potential for
various applications, including dynamic simulations, 3D printing, virtual
reality, and industrial design. By ensuring that the reconstructed objects
adhere to real-world physical laws, our method can enhance the realism and
functionality of virtual environments, improve the precision of 3D printed
objects, and contribute to the development of more reliable industrial
designs.

There are mainly two potential negative societal impacts: Improved 3D
reconstruction capabilities could potentially be misused to create highly
realistic fake objects or environments for disinformation purposes. This could
include generating deceptive media content that appears authentic. As the
framework automates the reconstruction process, there is a potential risk of
it being used in automated systems without sufficient oversight, potentially
leading to unintended and harmful outcomes due to errors or misuse. Developing
systems to monitor the use of the technology and ensure accountability for its
applications, as well as providing comprehensive guidelines and training for
users to promote ethical use and awareness of potential misuse, will address
these potential negative impacts.