  The key task of physical simulation is to solve partial differential
equations (PDEs) on discretized domains, which is known to be costly. In
particular, high-fidelity solutions are much more expensive than low-fidelity
ones. To reduce the cost, we consider novel Gaussian process (GP) models that
leverage simulation examples of different fidelities to predict
high-dimensional PDE solution outputs. Existing GP methods are either not
scalable to high-dimensional outputs or lack effective strategies to integrate
multi-fidelity examples. To address these issues, we propose Multi-Fidelity
High-Order Gaussian Process (MFHoGP) that can capture complex correlations both
between the outputs and between the fidelities to enhance solution estimation,
and scale to large numbers of outputs. Based on a novel nonlinear
coregionalization model, MFHoGP propagates bases throughout fidelities to fuse
information, and places a deep matrix GP prior over the basis weights to
capture the (nonlinear) relationships across the fidelities. To improve
inference efficiency and quality, we use bases decomposition to largely reduce
the model parameters, and layer-wise matrix Gaussian posteriors to capture the
posterior dependency and to simplify the computation. Our stochastic
variational learning algorithm successfully handles millions of outputs without
extra sparse approximations. We show the advantages of our method in several
typical applications.
