1. Introduction. The coronavirus disease 2019 (COVID- 19) and resulting pandemic has had a massive impact on worldwide public health (1.2 million deaths and 46 million reported infections as of 31 October 2020 1 ), nations' economies (a projected 4.5% reduction in global GDP in 2020 despite over $18 trillion in stimulus 2 ), and the downstream consequences on human welfare.
Since SARS-CoV-2, the virus underling COVID-19, is a novel virus, much remains uncertain about its transmission and spread. Computer modeling is thus fundamental to the projections needed to inform decision making and policy. A hallmark of the COVID-19 pandemic is that it is characterized by severe heterogeneities across populations: in infectiousness and susceptibility, in age-and co-morbidity-dependent morbidity and mortality rates, in human behavior and government policy, in medical care, and so on. Standard compartment-type population models in the form of systems of ordinary differential equations assume well-mixed homogeneous populations [29] and thus are poorly suited to modeling strongly heterogeneous disease spread.
Fully accounting for heterogeneities in infectious disease models, on the other hand, introduces numerous parameters. Solution of the resulting large scale inverse problems presents challenges due to the typically noisy, sparse, and inconsistently reported epidemiological data.
Moreover, due to the uncertainties in both models and data, quantification of uncertainties in parameter inference and subsequent predictions is paramount. Adopting the framework of Bayesian inference to account for these uncertainties leads to a high dimensional Bayesian inverse problem governed by (for example) differential equations. The twin curses of high dimensionality of parameter space and complexity of the forward model have historically made Bayesian inversion governed by differential equation models prohibitive. However, a number of theoretical and algorithmic advances over the past decade are making large-scale (partial) differential equation-based Bayesian inversion tractable, based on exploitation of the geometric structure and intrinsic low dimensionality of the posterior. These include geometrically-informed MCMC methods [3, 7, 20, 21, 22, 36, 38, 45] as well as more recent Stein variational methods [9, 10, 17, 23] . These methods have made Bayesian inversion feasible for complex problems in ice sheet dynamics and global seismology with as many as O(10 6 ) parameters [4, 6, 27] .
Our goal in this paper is to specialize such methods-in particular Stein variational methods, which are particularly appropriate for highly nonlinear parameter-toobservable maps-to high-dimensional heterogeneous infectious disease models, with specific application to modeling COVID-19 accounting for long-term care (LTC) facilities. A large percentage (over 40% 3 ) of COVID-19 deaths have occurred in nursing homes and other LTC facilities in the United States (and many other countries), despite the small percentage (∼0.7%) of the population in these facilities. Modeling, inference, and prediction under uncertainty of the spread and severity of COVID-19-while distinguishing between populations both inside and outside LTC facilities as well as their interactions-is crucial for implementing optimal strategies to contain and mitigate the disease through both non-pharmaceutical and pharmaceutical interventions. We endow our COVID-19 epidemiological model (see below) with a large number of uncertain parameters to account for variations in transmission, testing, hospitalization ratio, and fatality ratio over time, as well as other fixed parameters. This results in a challenging Bayesian inverse problem governed by ODEs with ∼1500 parameters, for which we employ the recently developed projected Stein variational gradient descent method [9] . The success of this method is predicated on the existence of an intrinsic low-dimensional subspace on which the parameters influence the observables. We investigate whether the inference problem for COVID-19 dynamics satisfies this property for our model.
There have been several different classes of epidemiological models developed to infer the dynamics of COVID-19 and make forecasts 4 , including classical compartmental models in the form of ordinary differential equations [34, 37, 41] and partial differential equations [28, 32, 43] , network models [18, 33] , statistical regression models [19, 42, 46] , agent-based models [24, 39] , etc. See [26, 29, 30] for general references on mathematical modeling of infectious diseases. Various strategies to account for heterogeneities in models of the dynamics of COVID-19 have been considered, including age, underlying health condition, connectivity, susceptibility, infectiousness, and spatial heterogeneity. However, given the significant expense of Bayesian inversion in high dimensions with complex models, a common approach to inference of parametrized models from data has been to invoke a low dimensional ansatz on the form of the parametric variation (such as assuming the variation in time of the basic reproduction number is represented by a sigmoid function), resulting in just a handful of parameters. One of our aims is to show that one can relax such assumptions to better capture the complex dynamics by accounting for the time dependent parameter uncertainties via stochastic processes, and attack the high dimensional inverse problem that results upon discretization with specialized methods such as projected Stein variational gradient descent.
Contributions. In this work, we propose to address the challenging ill posed inverse problem of Bayesian inference of COVID-19 models with high dimensional parameter spaces from noisy data, by exploiting the intrinsic low-dimensionality of the parameter-to-observable map and through the use of the recently developed projected Stein variational gradient descent method [9] , which requires a number of model solutions that is independent of the parameter dimension. We show that it is feasible to infer models of COVID-19 dynamics in very high parameter dimensions by aggressively exploiting problem structure. As a testbed application that is important in its own right, we model heterogeneity and severity of COVID-19 accounting for disease dynamics inside and outside LTC facilities and their interactions. We are not aware of any model in the literature that specifically address the heterogeneity associated with LTCs and the general population. We employ a heterogeneous compartmental model with the compartments of susceptible, exposed, infectious, recovered, hospitalized, deceased, confirmed, and unconfirmed cases to describe the epidemics both inside and outside LTC facilities. To better capture the dynamics, we consider a subset of model parameters to be continuously varying in time; these include (1) an effective virus transmission reduction factor that represents a compound effect of, e.g., social distancing, isolation, quarantining, and mask wearing; (2) an infection confirmation ratio that reflects the testing capacity and strategy; (3) an infection hospitalization ratio, and (4) a hospitalization fatality ratio, the latter two of which relate to time dependent changes in hospital capacity (ICU beds, ventilators), effectiveness of treatment, infection variations with age and health condition of the population, etc. This gives rise to an infinite dimensional Bayesian inverse problem [40] , which when discretized in time, results in a high dimensional problem.
To infer the model parameters, including the ∼1500 time-dependent parameters, we use reported data including the daily number of deceased cases from inside and outside LTC facilities separately, as well as the daily number of confirmed cases and the number of currently hospitalized cases for both inside and outside LTC facilities combined. We present both a deterministic inversion approach with a suitable regularization, and a Bayesian inference approach with proper prior distributions to learn these parameters with quantified uncertainties. We present the results of inference for both approaches, with data from both New Jersey and Texas that represent different development stages of COVID-19. Our model can capture well all reported data with 90% credible intervals covering the 7-day averaged data. Moreover, we present a future forecast and validation of the number of confirmed, infected, hospitalized, and deceased cases using our inference results. In a follow-on article, we use these posterior predictives as a basis for stochastic optimal control to determine optimal mitigation strategies for COVID-19 under uncertainty.
The paper is structured as follows. We present the heterogeneous compartment model in Section 2, followed by presentation of the available data and their post-processing in Section 3. Section 4 is devoted to presentation of a deterministic inversion formulation with appropriate regularization of the parameters and the inference results. This serves as a first step for a Bayesian inference formulation with proper prior distributions in Section 5. We introduce the projected Stein variational gradient descent method for the Bayesian inference problem and the results in Section 5. Finally, we draw some conclusions and perspectives in Section 6.
To learn the epidemics of COVID-19 both inside and outside longterm care facilities (LTC), we present a data-driven modeling of the epidemics, where the model is constructed by the principle of (1) being able to explain the available reported data, specifically the number of daily deaths inside and outside LTC for a certain period of time, the number of currently hospitalized cases, and the number daily confirmed cases, and (2) as simple as possible to avoid over-parametrization and overfitting. Under such principle, we consider the following compartmental model, as shown in Figure 1 , with two groups to describe the epidemics of COVID-19 inside and outside LTC within a time period t ∈ (t start , t end ]
where the compartmental variables S i , E i , I i , H i , R i , D i represent the number of susceptible, exposed, infectious, hospitalized, recovered, and deceased cases inside LTC (i = 1) and outside LTC (i = 2), and P c i and P u i represent the accumulated number of positive (infected) cases that are confirmed/reported by viral test or medical diagnosis, and unconfirmed/unreported inside LTC (i = 1) and outside LTC (i = 2). The effective transmission factor C i (t) is given by
Here the numbers of population inside LTC and outside LTC are represented by N 1 and N 2 . The time-dependent parameter α i (t) ∈ (0, 1) accounts for an effective reduction of the transmission at time t for group i due to, e.g., social distancing, isolation, and quarantine. The contact matrix C ∈ R 2×2 represents the number of contacts inside (C 11 ) and outside (C 22 ) LTC, as well as the interactions between them (C 12 and C 21 ) through LTC and hospital staff, and family visitors to LTC. Specifically, C 12 represents the number of contacts between susceptible LTC residents and infectious others, while C 21 represents the number of contacts between infectious LTC residents and susceptible others. The parameters β i , σ i , η i , µ i , γ I i , γ H i are the transmission rate, the latency rate, the hospitalization rate, the deceased rate, and the recovery rate from infectious and hospitalized stages, which are the inverse of the transmission period, the latency period (exposed but not infectious yet), from infectious to hospitalized period, from hospitalized to deceased period, from infectious to recovery (not infectious) period, and from hospitalized to recovery (leaving hospital) period. The parameter π i (t) is related to the infection hospitalization ratio (IHR) denoted as ζ i (t) ∈ (0, 1), i.e., the proportion of the infected cases that are hospitalized
where we assume that the IHR is time dependent possibly due to (1) the dynamic change of the severity of the infected population (e.g., age and health condition) and (2) hospital availability. Similarly, the parameter ν i (t) is related to the hospitalization fatality ratio (HFR) denoted as ξ i (t) ∈ (0, 1), i.e., the proportion of the hospitalized cases that are deceased
where we assume that the HFR is also time dependent possibly due to (1) the dynamic change of the severity of the hospitalized population (e.g., age and health condition) and (2) medical resources (e.g., ICU bed, ventilator, and new drugs). Finally, the parameter τ i (t) represent the infection confirmation ratio (ICR), or the proportion of infected cases that are confirmed/reported positive by viral test or medical diagnosis, which depends on the test strategy and capacity. 31 . The LTC share of New Jersey is close to the US average and that in many European countries. The LTC share of Texas is close to 50% early on and decreasing to 31% when the LTC became strongly protected. Moreover, the two states have distinctive development of the epidemics: New Jersey experienced fast outbreak in an early stage (March and April) with a large number of deaths, and has since flatten the curve by implementing strict mitigation policy of lockdown, social distancing, etc., while Texas had slower spread of the virus than New Jersey early on (March and April) with a much smaller number of deaths, yet turned to an accelerated increase of cases after quickly relaxing the mitigation (from May to July) and flattened the curve since mid-July. These epidemic developments can be observed in Figure 2 for New Jersey and Figure 3 for Texas. The number of COVID deaths inside and outside LTC by state of US is reported in Nursing Home/LTC COVID Deaths By State 5 for a certain period of time. The bottom-right part of Figure 2 displays the reported number of accumulated and daily deaths inside LTC in New Jersey from April 17 to August 31. The number of daily deceased cases has some negative values and is very oscillatory showing a 7 day periodic pattern with much smaller number of deaths reported during the weekend than during the week. We treat this oscillatory pattern as report error due to report delay, and possibly incompleteness and incorrectness. To correct the report error and make the daily numbers smoother in a given time interval (t 1 , t 2 ), we compute a seven day moving average d r i , i = 1, 2, from the reported daily datad r i , as
The smoothed data by seven day moving average are also shown in Figure 2 . Note that the number of daily deceased cases is already relatively big on the first day of report, which can not reflect earlier stage of the epidemics. To learn the early dynamics, we use another set of data from the COVID19 Tracking project 6 , as shown in the other three figures in Figure 2 , which reports the number of daily deceased cases, infected cases that are confirmed by viral test or medical diagnosis from early March, as well as the number of currently hospitalized cases. Note that these data are aggregated for two groups of inside and outside LTC, which are not available separately. To smooth these oscillatory data, we also compute a seven day moving average p r (t), d r (t), and H r (t) for the number of daily confirmed casesp r (t), daily deceased casesd r (t), and currently hospitalized casesH r (t). As the two sets of data are slightly inconsistent for the reported number of deceased cases, we use d r 1 (t) for the number of deaths inside LTC, and d r 2 (t) = d r (t) − d r 1 (t) as that outside LTC to keep a consistent d r (t) in the entire time period. We denote the number of accumulated confirmed and deceased cases for the smoothed data as 
As the compartmental variables of the dynamic model (2.1) is very sensitive to the parameters due to their exponential growth in time, the uncertainty on the (pre-inferred) parameters would lead to large uncertainty or even blowup in the model output. To avoid such blowup, we first deterministically infer the parameters at which the model can optimally fit the data in this section, and then quantify the uncertainty for the inferred parameters in next section. We formulate the deterministic inversion problem, and present the results for New Jersey and Texas.
To infer the parameters in (2.1) given data presented above, we minimize an objective function with respect to the parameters, which consists of a misfit function between the model output and the data, and a penalty function that imposes certain constraints on the parameters.
By D(t) at t ∈ [t start , t end ] we denote the number of total deaths from the model (2), given as
Then we compute the number of daily deaths from the model (2) at t ∈ [t start +1, t end ] as
Similarly, we denote P c (t) as the number of total infected cases that are confirmed, Figure 2 for Texas here. The increase of the number of confirmed, deceased, and hospitalized cases in Texas was slower than that of New Jersey in the early stage, and got accelerated before flattened lately.
given by
and the number of daily confirmed cases as
As the early data (deaths and confirmed cases) are subject to randomness of the epidemics in the early stage, we use the data from the time when the numbers are relatively big. Let t p 1 denote the first day that the number of total confirmed cases surpassed 100, t h 1 denote the first day that the number of currently hospitalized cases surpassed 10, t d 1 denote the first day that the number of total deaths surpassed 10, t d 2 denote the day before the number of LTC deaths is reported, and t end the last day the data is available. To this end, we define the misfit function as
where we use a small weight 0.1 for the daily data as they are more oscillatory thus hard to match and probably less reliable than the accumulated data. Moreover, these weights are roughly proportional to the inverse of the standard deviation of the data (the difference between the reported number and seven day averaged number) that reflects the uncertainty or reliability of the reported data. The aggregated parameter set θ in (4.5) collects all the parameters appearing in model (2.1) as (4.6)
where t set = (t 0 , t 1 , . . . , t k ) is a sequence of equispaced time points with step size δt.
Such an interpolatory definition with a relatively big step size, e.g., δt = 7, is employed to prevent over-parametrization that may lead to daily oscillatory and large variation in the transmission reduction effect, which plays a role of regularization in the deterministic inversion by solving a relatively low-dimensional optimization problem. The optimizer can be interpolated and used as an initial guess to refine the optimization with smaller step size, e.g., δt = 1, which is often easier than directly solving a high-dimensional optimization problem. In the misfit function defined in (4.5), we use data of the number of both accumulated and daily confirmed and deceased cases, as we find that the accumulated data can prevent over-fitting for the oscillatory (even after smoothing by moving average) daily data. Since these numbers are different by several orders of magnitude, we use their logarithmic scale to make them comparable. Moreover, we explicitly write θ as one argument for the quantities in the misfit function (4.5) that depends on θ through the model (2.1).
We consider the following penalty for the parameters
where the first term in the first line represents the change of the transmission reduction factor α i (t) in time t, an approximation to the variational integral t k t0 (dα i (t)) 2 dt by a trapezoidal rule, which is penalized to prevent its large variation; the second term in the first line is used to penalize decreasing α 1 (t), which reflects an increasing transmission reduction inside LTC as strict mitigation policy has been implemented and likely maintained in preventing excessive death that has accounted for a significant proportion of total deaths despite a small population, about 0.7% of total population in average in US. Similar penalizations on the change of the ICR, IHR, HFR in time are considered. In addition, we consider a penalization on the deviation of ICR, IHR, HFR from their reference values weighted by a scaling parameter s t > 0, which is set to be large that reflects less confidence and allows relatively large deviation from their reference values. In the last two lines of of the penalty function (4.7), the devi-
Given the misfit function (4.5) and the penalty function (4.7), we define an objective function
where λ > 0 is a scaling parameter that balances the two terms. The deterministic inversion problem can then be formulated as the optimization problem: find θ * ∈ Θ, such that
where Θ is an admissible parameter range. We start the optimization with a relatively big step size δt = 7, which leads to a 224-dimensional parameter θ for both New Jersey and Texas, including dimension 26 for each of the time dependent variables, and dimension 16 for the rest of the time independent parameters. To refine the optimization with smaller step size δt = 1, we use the interpolated optimizer as the initial guess, which leads to a 1520-dimensional parameter θ for New Jersey (1488 for Texas), including dimension 188 for New Jersey (184 for Texas) for each of time dependent variables. To solve the high-dimensional optimization problem (4.9) with bounded constraint on the parameter θ ∈ Θ, we use a Python implementation of an optimizer SLSQP, a Sequential Least SQuares Programming algorithm, which is a gradient-based quasi Newton method with a BFGS update for the approximation of the Hessian [31] . The gradient of the objective function J(θ) with respect to the high-dimensional parameter θ is computed by automatic differentiation using an Autograd 7 library. We also develop an accelerated computation of the gradient by an adjoint approach, for which the expressions of involved derivatives are first symbolically computed by a symbolic library SymPy 8 and then accelerated by a high performance Python compiler Numba 9 . See an abstract formulation of the adjointbased computation of the gradient in Appendix A, which is used in next section for Bayesian inference, about 100X faster than the automatic differentiation by Autograd.
. We present the setup and the results of the deterministic inversion for New Jersey and Texas. We initialize the model (2.1) with the number of susceptible cases S 1 and S 2 the same as the population inside and outside LTC; we set E 1 = 1 exposed case inside LTC and E 2 = 100 exposed cases outside LTC, and all other compartments as zero. To initialize the parameters, we set α 1 (t) = α 2 (t) = 0, , t ∈ t set , and set a bound as α i (t 0 ) ∈ [0, 0.1] to have a small transmission reduction at the beginning, and α i (t) ∈ [0, 0.9], t = t 1 , . . . , t k , where the upper bound 0.9 assumes that the transmission can not be completely prevented by the mitigation policy. For the reference values of the parameters presented in the following, besides a few assumptions, we consider those available from [24, 33, 42] , CDC COVID-19 Pandemic Planning Scenarios 10 and references therein. More specifically, we initialize ICR τ 1 (t) = 40% and τ 2 (t) = 10%, IHR ζ 1 (t) = 25% and ζ 2 (t) = 5%, HFR ξ 1 (t) = 40% and ξ 2 (t) = 10%, which lead to an initial infection fatality ratios (IFR = IHR×HFR) 10% inside LTC and 0.5% outside LTC. We also use these values as the reference values in the penalty term (4.7). We consider the bound [ρ/4, 2ρ] ⊂ [0, 1], withρ representing the initial values of the ratios τ i (t), ζ i (t), ξ i (t). For the rest of parameters in θ, we initialize them at their reference valuesρ as reported in Table 1 , and consider the bound [ρ/2, 2ρ] to reflect relatively large uncertainty. Specifically, for the latency rate we use σ 1 = σ 2 = 1/4 (inverse of the latency period, exposed but not infectious yet, 4 days), for hospitalization rates we use η 1 = η 2 = 1/8 (inverse of the period from hospitalization to death, 8 days), for the recovery rates we use γ I 1 = γ I 2 = 1/5 (inverse of the infectious period, 5 days) and γ H 1 = γ H 2 = 1/7 (inverse of the hospital stay period, 7 days). For the contact matrix we use C 11 = 4, C 12 = 0.5, C 21 = 0.1, C 22 = 2, which reflect relatively more frequent contacts inside LTC than outside LTC, and larger infection influence of people from outside to inside LTC than those from inside to outside LTC. The contact matrix is rescaled by the transmission parameters β 1 = β 2 = 0.3 such that the initial doubling time of infection and death is about 2 ∼ 3 days.
We solve the optimization problem (4.9) for deterministic inversion of the parameters with the scaling parameter λ = 100 to balance the misfit and penalty terms. Table 1 Reference values and optimal values of the time independent parameters in θ collected in (4.6).
The optimal values are obtained by solving the deterministic inversion problem (4.9) θ reference values optimal values (NJ) optimal values (TX) # daily deceased data model Fig. 4 . The same as in Figure 2 . The data here is the same as the seven day moving average data in Figure 2 . The model output is the result of model (2.1) at the optimal parameter value obtained by deterministic inversion.
We use relatively large scaling parameters s t = 10, s = 5 in the penalty term (4.7) such that the optimal parameters can be easily deviated from the reference values. We run the SLSQP optimization algorithm and stop it at 100 iterations, when the objective function (4.8) does not decrease much (less than 10 −4 compared to the decrease at the first step). We plot the comparison between the model output and the smoothed data (with seven day moving average) in Figure 4 for New Jersey and Figure 5 Figure 3 . The data here is the same as the seven day moving average data in Figure 3 . The model output is the result of model (2.1) at the optimal parameter value obtained by deterministic inversion.
be observed for both states, thanks to the flexibility of the parametrization especially in the time dependent transmission reduction factors, infection confirmation ratios, infection hospitalization ratios, and hospitalization fatality ratios for both inside and outside of LTC, whose optimal values are shown in Figure 6 for New Jersey and Figure 7 for Texas. This demonstrates the effectivity of both the parametrization and the optimization algorithm in the framework of deterministic inversion. We remark that some data and model mismatch can be observed, e.g., the LTC deceased cases in both states can be probably accounted by the poor quality of the data due to delay or misreport, which can be supported by the very oscillatory (already smoothed by 7 day moving average) and the abruptly changed (in early July) number of daily deaths especially in Texas. This slight mismatch can also be explained by the use of the relatively small weight for the daily data in the misfit function (4.5), which is to prevent catching the artifact of oscillation and overfitting that may lead to too much and unrealistic oscillation in the time dependent parameters. Nevertheless, the data and model output are matched rather well for the number of accumulated confirmed and deceased cases. We also remark that the data and model match for the number of hospitalized cases is very high, even for the sharp change, e.g., in mid-April and late August in New Jersey, and in early April, mid-June and late July in Texas.
We make a few observations and comments for the deterministic inversion results of parameters. The optimal transmission reduction factor inside and outside LTC are shown in the top-left part of Figure 6 for New Jersey and Figure 7 for Texas. We can see a sharp increase in late March and early April for both inside and outside LTC in both states, when increasingly stricter mitigation policy was implemented, e.g., school closure, prohibitive large group gathering, lockdown, limited or prohibitive family visiting to LTC, increasing use of personal protective equipment (PPE) inside LTC. The transmission reduction is relaxed by the gradual lift of the strict mitigation measures with several phases starting from May. However, it is observed that the lift was faster in Texas than in New Jersey in May and June outside LTC, which led to the sharp increase of the number of confirmed and hospitalized cases in Texas, and continuing decrease in New Jersey. The transmission reduction factor is increased in July in Texas probably because the lift of effective mitigation measures were reversed, while wearing masks was suggested and later mandated. We remark that the relatively more effective transmission reduction inside LTC than outside has been maintained probably by large-scale testing and widely use of PPE. We observe that the infection hospitalization ratios and hospitalization fatality ratios outside LTC in both states has decreased from their initial values, which may be explained by the shift of the new infections from the general population to the younger and healthier population, which is built into the penalization (4.7). We remark that the unexpected increase of the hospitalization fatality ratio in Texas from mid-July is probably due to the change of the reported fatality from July 27 and on. 11 The optimal values for the other time independent parameters are reported in Table 1 , some of which are evidently different from the reference values. The slightly different values between New Jersey and Texas imply slightly different transmission/confirmation/hospitalization/death/recovery pattern. However, we caution that these optimal values may differ from the reality even though they provide the best fit for the data because (1) the inference problem is ill-posedness, i.e., the data are noisy, limited, and can only effectively inform the most data sensitive parameters, (2) the model is inadequate that may not capture the heterogeneity in age and health condition, different susceptibility, infectivity, and symptoms of different people, and (3) the partial information embedded in the penalization on the model parameters may not represent the reality.
In learning the epidemics through assimilating the available data to the model (2.1) by solving the optimization problem (4.9), various inevitable uncertainties are neglected, including the parameter uncertainty and data uncertainty, besides the uncertainty in the model itself. Without accounting for these uncertainties, the parameters inferred from the smoothed data may not capture the reality, which may lead to incorrect forecast results for the future transmission of the virus and ineffective optimal mitigation design based on the deterministic inversion. In this section, we deal with the parameter and data uncertainties in a Bayesian framework for statistical inference of the parameters equipped with probability distributions and provide confidence for the model output. In particular, to solve the Bayesian inference problem with high-dimensional and heterogeneous parameters (time-dependent and time-independent), we employ our recently developed pSVGD method [9] , which exploits the intrinsic low dimensionality of the parameter-to-observable map and efficiently draw samples from the posterior distribution by solving an optimal transport problem with parameter projection into low-dimensional subspaces.
5.1. Formulation. The parameter θ defined in (4.6) is a set of heterogeneous parameters, including the discrete sequences of the four types of time dependent parameters (α i (t)) t∈tset , (τ i (t)) t∈tset , (ζ i (t)) t∈tset , (ξ i (t)) t∈tset , i = 1, 2, the time independent parameters β i , σ i , η i , µ i , γ I i , γ H i , i = 1, 2, and a contact matrix C. We present prior probability distributions for these parameters. For simplicity, we assume that these parameters are independently distributed.
We construct the prior distributions for (τ i (t)) t∈tset , (ζ i (t)) t∈tset , (ξ i (t)) t∈tset in the same way as for (α i (t)) t∈tset , and only present it for (α i (t)) t∈tset for simplicity. The sequences (α i (t)) t∈tset are obtained as a temporal discretization of the transmission reduction factor α i (t) for t ∈ (t 0 , t k ). To equip a prior probability distribution for α i (t) ∈ (0, 1), we assume that it is given by the form
where tanh is the hyperbolic tangent function so that the constraint α i (t) ∈ (0, 1) is always satisfied. Moreover, the time-dependent function g i (t) ∈ N (g * i , C g ) is assumed to be a Gaussian process with mean g * i and covariance C g . We set the mean as
is defined as piecewise linear interpolation of the optimal sequences (α * i (t)) t∈tset , which is obtained as the solution of the optimization problem (4.9) for deterministic inversion. We assume the covariance C g is given by
where B is an elliptic differential operator in time domain (t 0 , t k ) given by
equipped with homogeneous boundary condition, where s g > 0, s I > 0 are scaling parameters, ∆ t is the Laplacian in time, and I is identity. Such a covariance (5.3) is often used in defining Gaussian process/random field with s g and s I controlling its correlation length and variance. By a temporal discretization at equispaced points t 0 , . . . , t k using finite difference, we have
for any twicely differentiable function g : [t 0 , t k ] → R, we obtain a discrete Laplacian L with homogeneous Neumann boundary condition, i.e., g (t 0 ) = g (t k ) = 0, such that
for g = (g(t 0 ), . . . , g(t k )) T , where the first term has the same form as that in the penalty term (4.7) up to a constant. To this end, we obtain the prior distribution for the vector g i = (g i (t 0 ), . . . , g i (t k )) T from the discretization of the Gaussian process g i ∼ N (g * i , C g ), as g i ∼ N (g * i , C g ) with the mean vector g * i = (g * i (t 0 ), . . . , g * i (t k )) and covariance matrix C g = (−s g L + s I I) −1 . By slight abuse of notation, we denote g as a collection of all the 4 × 2 (4 for different ratios, 2 for inside and outside LTC) time dependent parameters, with a total dimension of 8(k + 1). Meanwhile, by slight abuse of notation we denote C g as a block diagonal covariance matrix, where each of the eight blocks is the covariance matrix (−s g L + s I I) −1 .
For the time independent parameters, including the scalar parameters of various rates β i , σ i , η i , µ i , γ I i , γ H i , i = 1, 2, and the contact matrix C, we collect them as a vector of parameters (5.7) f = (β 1 , β 2 , σ 1 , σ 2 , η 1 , η 2 , µ 1 , µ 2 , γ I 1 , γ I 2 , γ H 1 , γ H 2 , C 11 , C 12 , C 21 , C 22 ) T , and assume that it has a prior distribution as
where h * = log(f * ) for f * defined as the optimal values of f obtained from the deterministic inversion problem (4.9). For the covariance C h , we assume that it is diagonal with the diagonal entries given as the entries of the vector (s h h * ) 2 with a scaling parameter s h > 0. Under this assumption, we have
which has the same form as that in the penalty term (4.7) up to a constant. Finally, we collect all the parameters as x = (g, h) ∈ R dx of dimension d x . Then the prior density of x, denoted as p 0 (x), satisfies
Note that − log(p 0 (x)) has similar form as part of that in the penalty term (4.7). However, the parameter x here, after the nonlinear transformation by the function in (5.1) and the logarithmic function in (5.8), is different from θ in (4.7). We remark that the deterministic inversion problem (4.9) can also be formulated as an optimization problem with respect to the transformed parameter x. However, we observe from tests that it is much more difficult to solve than the one with respect to the parameter θ, mainly due to the stiffness of the nonlinear transformation, especially (5.1), in a major part of the whole parameter space. Nevertheless, by defining the prior distribution of x centered at the transformed optimal values of the parameter θ, and use a proper covariance for x, the stiffness can be effectively alleviated in the Bayesian inference problem. Based on the above definition of the transformed parameter x and its prior distribution, we assume that the data and model output can be formulated as
where the data y represents a collection of all the data as given in the misfit function and the parameter-to-observable map u(x) is given by
where each set of data in y and model output in u(x) has a time span corresponding to that in (4.5) . Note that in u(x) the model output depends on the parameter x through θ by the transformation (5.1) and (5.8) , as well as the model (2.1) which depends on the parameter θ. For simplicity, we assume that the additive noise ω ∼ N (0, C ω ) with zero mean and a symmetric positive definite covariance matrix C ω . Then the likelihood function, denoted as f (x), satisfies
By Bayes' rule, the posterior density of the parameter x is given by
where Z is a normalization constant given by
which is typically intractable to compute, especially for a large dimension n. The tasks of Bayesian inference is to draw samples of the parameter x from its posterior distribution, and compute some statistical quantity of interest, e.g., mean, variance, and confidence interval of the parameter and the model output. Several challenges are faced for these tasks, including (1) curse of dimensionality -the computational complexity grows rapidly (typically exponentially) for conventional methods with respect to the parameter dimension, and (2) the geometry of the posterior density is complex, e.g., non Gaussian, multimodal, concentrating in local regions, etc. We have developed several approaches to address the curse of dimensionality, including [9, 10, 12, 13, 14, 15, 17] that exploit the sparsity, low dimensionality, and geometry of the parameter-to-observable map. In this work, we employ our recently developed projected Stein variational gradient descent (pSVGD) method [9] , which can break or considerably alleviate the curse of dimensionality by exploiting the intrinsic low-dimensionality of the inference problem using gradient information matrix, and effectively deal with the complex geometry by transporting samples in parallel from the prior distribution to the posterior through optimization.
We briefly present SVGD from [35] and then introduce the projected SVGD from [9] to tackle high-dimensional inference problems. As a variational inference method, SVGD seeks an approximate density q ∈ Q in a density function space Q for the posterior density p by minimizing a Kullback-Leibler (KL) divergence that measures the difference between two densities, i.e., (5.17) q
where D KL (q|p) = E x∼q [log(q/p)], which vanishes when q = p. A specific function class is Q = {T p 0 : T ∈ T }, where T is a pushforward map that pushes the prior density to a new density q := T p 0 through a transport map T (·) : R dx → R dx defined in a space T . For
where φ : R dx → R dx is a differentiable perturbation map, and > 0 is a learning rate, there holds
where A p is the Stein operator given by
A practical SVGD algorithm was obtained by choosing T = (H dx ) dx = H dx ×· · ·×H dx , a tensor product of a reproducing kernel Hilbert space (RKHS) H dx with kernel k(·, ·) : R dx × R dx → R, which push the samples x 0 1 , . . . , x 0 N drawn from the prior p 0 step by step as (5.21) x +1 m = x m + lφ * (x m ), m = 1, . . . , N, = 0, 1, . . .
whereφ * (x m ) is an approximate steepest direction given by
The kernel k plays a critical role in pushing the samples to the posterior. One choice is Gaussian
where h x is a bandwidth, e.g., h = med 2 / log(N ) with med representing the median of sample distances [35] . However, it is known that the kernel suffers from the curse of dimensionality for large dimension d x , which leads to the samples that are not representative of the posterior, as observed in [9, 17, 44, 47] .
To address the curse of dimensionality, we exploit the property that the subspace in which the posterior differs from the prior is low dimensional due to the ill-posedness of the inference problem, as proven for some model problems and numerically demonstrated for many others [1, 2, 5, 8, 9, 15, 16, 17, 25, 27, 36] . We project the parameter x into this subspace spanned by r basis functions X r = span(ψ 1 , . . . , ψ r ) as
where the matrix Ψ r := (ψ 1 , . . . , ψ r ) ∈ R dx×r represents the projection matrix and the vector w := (w 1 , . . . , w r ) T ∈ R r is the coefficient vector with element w i := ψ T i x for i = 1, . . . , r. The basis functions ψ 1 , . . . , ψ r are obtained as the eigenvectors for the generalized eigenvalue problem (5.25) Hψ
represent the dominant eigenpairs of (H, C x ) correspond to the r largest eigenvalues λ 1 ≥ · · · ≥ λ r . The matrix C x is the prior covariance of the parameter x = (g, h) , i.e., C x = diag(C g , C h ), a block diagonal matrix with each diagonal block representing the covariance of g and h. The matrix H is chosen as a gradient information matrix, which is defined as the average of the outer product of the gradient of the log-likelihood w.r.t. the posterior, i.e.,
which is evaluated by sample average approximation to be specified later. By the projection (5.24), we consider a decomposition of the prior for the parameter x = x r + x ⊥ , where x r = P r x, as
where p r 0 and p ⊥ 0 are marginal and conditional densities for x r and x ⊥ . Moreover, we define a marginal likelihood function (an optimal approximation of the full likelihood function [9] ) as
where X ⊥ is the complement of the subspace X r . Since the uncertainty in x r = P r x = Ψ r w can be fully represented by the coefficient w, we consider the prior and posterior density for w, given by
where the normalization constant Z w = E w∼π0 [g(Ψ r w)].
To sample from the posterior π in (5.29), we employ the SVGD method presented above for the projection coefficient. Specifically, pSVGD pushes the samples w 0 1 , . . . , w 0 N from π 0 (w) as ∇ w n log π(w n )k r (w n , w m ) + ∇ w n k r (w n , w m ),
where the kernel k r is specified as
where the metric matrix Λ+I with eigenvalues Λ = diag(λ 1 , . . . , λ r ) of (5.25) are used to account for data impact in different directions ψ 1 , . . . , ψ r . The pSVGD algorithm is presented in Algorithm 5.1, which is implemented in parallel with samples distributed in multiple cores. We remark that the projected likelihood in (5.29) can be approximated by its definition (5.28) as
set the observation noise covariance matrix C ω as the identity matrix. This leads to sufficient variation of the parameters that result in large variation (and no blowup) of the model output at different samples drawn from the prior distribution, as can be observed in the left part of Figure 12 and 14 for New Jersey and Texas. We run the adaptive pSVGD Algorithm 5.2 with N = 125 samples in each of K = 8 cores, so 1000 samples in total. We stop the algorithm at 10 × 10 iterations, where the projection basis functions Ψ x r are updated for each 10 iterations. We remark that the move of the samples already becomes small towards the end of the iteration. We choose the dimension r such that the eigenvalue λ r+1 < 10 −1 . The decay of the eigenvalues at different pSVGD iterations are shown in Figure 8 for both states, from which we can observe a significant reduction of the eigenvalues by several orders of magnitude for the first few dimensions, which implies that the data informed parameter subspace is intrinsically low-dimensional, e.g., the first 10 dimension (out of the total 1520 for New Jersey and 1488 for Texas) catches over 99% difference of the prior and posterior distribution measured by the KL divergence (5.17) . Figure 9 displays the statistics of the posterior samples, including their mean and 90% credible interval, in which 90% samples fall, as well as the optimal values obtained from the deterministic inversion for the time independent parameters, while Figures 10 and 11 display these quantities for the time dependent parameters. We observe that the posterior sample mean is close to the optimal values obtained by the deterministic inversion. Meanwhile, it is evident from the large 90% credible interval of the posterior samples that there are uncertainties in the parameters that can not be ignored. These uncertainties become crucial for reliable forecast of new infected, hospitalized, and deceased cases into the future based on the model output, as well as for risk-averse optimal mitigation design to flatten the curve.
Finally, we show the comparison of the statistics (mean and 90% credible interval) of the model outputs at the prior samples and the posterior samples in Figures 12 and  14 for confirmed and hospitalized cases, and in Figures 13 and 15 for deceased cases. From these figures we can observe that the variation of the model outputs is considerably reduced from the prior to the posterior samples. Moreover, most of the data, except for a few data points that is largely uncertain or possibly incorrect/outlier, fall inside the 90% credible interval of the model outputs, which demonstrate the reliability and advantage of the Bayesian inference that account for the parameter hospitalized, and daily confirmed and deceased cases, while the credible intervals for the number of LTC deceased cases are relatively larger than that of the total number of deceased cases, etc., which may bring in over or under confidence on the inference results. This is probably due to the relatively simple choice of the prior distribution of the parameters and the distribution of the observation noise.
In this section, we compare the model forecast under uncertainty and compare the forecast results with the data. We divide the data into two parts, one part including all but the last four weeks, which are used for parameter inference, and the other part including only the last four weeks, which are used for comparison with the forecast. We run pSVGD for Bayesian inference of the parameters with 1000 samples for 100 iterations. We use the 1000 samples at the end of iteration to make forecast for the last four weeks, during which the time dependent parameters are assigned with values at the end of the inference period, i.e., one day before the last four weeks. The inference and forecast results for the confirmed, hospitalized, and deceased cases are shown in Figure 16 for New Jersey and Figure  17 for Texas. We can observe a relatively good agreement of most of the data and the model forecast at both the optimal parameters by deterministic inversion and the sample posterior mean. However, some of the model forecast evidently depart from the data, e.g., the hospitalized cases in Texas with sharp change, due to the evident dynamic change of epidemics reflected by the time dependent parameters (e.g., IHR and HFR outside LTC) in the last four weeks, see Figure 11 . Nevertheless, the data still fall in the 90% credible interval of the forecast. We can also observe that the longer the forecast is, the larger the uncertainty of the forecast results become.
6. Conclusions. We have proposed a Bayesian inference framework for learning the dynamics of COVID-19, with application to the spread and outcome inside and outside LTC facilities, which have experienced about 40% of total COVID-19 caused and/or related deaths. We designed an LTC stratified SEIR compartmental model with extended compartments of confirmed and unconfirmed positive infections, hospitalization, and death to account for the available data on the number of positively confirmed, hospitalized, and deceased cases, as well as the number of deceased cases inside LTC facilities. To estimate the parameters in the model, including the timedependent transmission reduction factors, various epidemic rates and proportions, as well as the contact matrix-and more importantly quantify the uncertainties in the parameters, model outputs, and data-we investigated both a deterministic inversion and a Bayesian inference approach that builds on the deterministic result. The deterministic inversion with our model was demonstrated to effectively capture the main dynamics, while the Bayesian inference method was demonstrated to achieve reliable quantification of various uncertainties, in particular providing tight 90% credible intervals that sufficiently cover the data. The high dimensional Bayesian inversion was made tractable by our pSVGD method's ability to exploit the two orders of magnitude effective dimensionality reduction inherent in inferring the model parameters Fig. 13 . The same as in Figure 12 for total deceased cases (top) and deceased cases inside LTC (bottom) in New Jersey.
from the available data. The posterior predictives produced by the Bayesian inference play a critical role in a solution of a risk-averse stochastic optimal control problem to find optimal mitigation strategies for COVID-19, which are presented in a follow-on article [11] .
Limitations of this work include: (1) The data we use inevitably include incorrect and delayed numbers, as discussed for the data in Figure 2 and 3. We accounted for these flaws only by smoothing the data using seven day moving averages. If they were available, additional data could improve the results of the inferences, e.g., the number of daily confirmed cases inside LTC, state dependent hospitalization ratios, time-dependent positive/negative ratio of the viral test, etc. (2) The model we constructed is, as mentioned in Section 2, as simple as possible for our target of learning the dynamics inside and outside LTC. Though it would have been straightforward to incorporate in our model, we did not account for the age-and risk-dependent heterogeneities in infection hospitalization and fatality ratios; their effect can be observed in the widening misfit in the number of deaths in recent months, probably due to the shift of infections toward a younger and healthier population. Figure 16 for Texas here.