Parallel numerical tensor methods for high-dimensional PDEs
A. M. P. Boelens D. Venturi D. M. Tartakovsky








stanford]A. M. P. Boelens
boelens@stanford.edu
ucsc]D. Venturi
venturi@ucsc.edu
stanford]D. M. Tartakovsky
tartakovsky@stanford.edu


[stanford]Department of Energy Resources Engineering, Stanford University, Palo Alto, CA 94305
[ucsc]Department of Applied Mathematics and Statistics, UC Santa Cruz, Santa Cruz, CA 95064


High-dimensional partial-differential equations (PDEs) arise in a 
number of fields of science and engineering, where they are used to describe the evolution of joint probability functions. Their examples include the Boltzmann and  Fokker-Planck 
equations. We develop new parallel algorithms to solve high-dimensional PDEs.
The algorithms are based on canonical and hierarchical numerical 
tensor methods combined with alternating least squares
and hierarchical singular value decomposition. Both implicit 
and explicit integration schemes are presented 
and discussed. We demonstrate the accuracy and 
efficiency of the proposed new algorithms in computing 
the numerical solution to both an advection 
equation in six variables plus time and a linearized version of the 
Boltzmann equation.



[2010] 00-0199-00







§ INTRODUCTION


High-dimensional partial-differential equations (PDEs) arise in a 
number of fields of science and engineering. 
For example, they play an important role in modeling 
rarefied gas dynamics <cit.>, 
stochastic dynamical systems <cit.>, 
structural dynamics <cit.>,  
turbulence <cit.>,  
biological networks <cit.>, and 
quantum systems <cit.>.

In the context of kinetic theory and stochastic dynamics, 
high-dimensional PDEs typically describe the evolution of 
a (joint) probability density function (PDF) of system states, 
providing macroscopic, continuum descriptions of Langevin-type stochastic 
systems. Classical examples of such PDEs include the Boltzmann 
equation <cit.> and the Fokker-Planck 
equation <cit.>. More recently, PDF equations 
have been used to quantify uncertainty in model predictions 
<cit.>.  

High dimensionality and resulting computational complexity of 
PDF equations can be mitigated by using particle-based methods 
<cit.>. 
Well known examples are direct simulation Monte-Carlo (DSMC) 
<cit.> and the Nambu-Babovsky 
method <cit.>. Such methods preserve physical 
properties of a system and exhibit high computational efficiency 
(they scale linearly with the number of particles), in particular in 
simulations far from statistical equilibrium <cit.>. 
Moreover, these methods have relatively low memory requirements 
since the particles tend to concentrate where the distribution function 
is not small. However, the accuracy of particle methods may 
be poor and their predictions are subject to significant statistical fluctuations 
<cit.>. Such fluctuations need to be post-processed  appropriately, for 
example by using variance reduction techniques. Also, 
relevance and applicability of these numerical strategies to 
PDEs other than kinetic equations are not clear.

To overcome these difficulties, several general-purpose algorithms 
have been recently proposed to compute the numerical solution 
to rather general high-dimensional linear PDEs. The most efficient techniques, 
however, are problem specific <cit.>. 
For example, in the context of kinetic methods, there is an extensive literature 
concerning the (six-dimensional) Boltzmann equation and its solutions close to 
statistical equilibrium. Deterministic methods for solving such problems 
include semi-Lagrangian schemes and discrete velocity models. 
The former employ a fixed computational grid, account for transport features of the 
Boltzmann equation in a fully Lagrangian framework, and usually adapt operator 
splitting. The latter employ a regular grid in velocity and a discrete collision
operator on the points of the grid that preserves the main physical properties
(see 3 and 4 in Di Marco and Pareschi <cit.> for an in-depth
review of semi-Lagrangian methods and discrete velocity models for kinetic
transport equations, respectively). 

We develop new parallel algorithms 
to solve high-dimensional partial differential equations
based on numerical tensor methods <cit.>. 
The algorithms we propose are based on canonical and 
hierarchical tensor expansions combined 
with alternating least squares and hierarchical 
singular value decomposition. The key element that opens 
the possibility to solve high-dimensional PDEs numerically with 
tensor methods is that tensor approximations proved to be 
capable of representing function-related N-dimensional data arrays 
of size Q^N with log-volume complexity 𝒪(Nlog(Q)). 
Combined with traditional deterministic numerical 
schemes (e.g., spectral collocation method <cit.>, 
these novel representations allow one to compute the 
solution to high-dimensional PDEs using low-parametric 
rank-structured tensor formats.

This paper is organized as follows. In Section <ref> we
discuss tensor representations of N-dimensional functions. In particular, 
we discuss canonical tensor decomposition and hierarchical tensor networks, 
including hierarchical Tucker and tensor train expansions. We also address the 
problem of computing tensor expansion via alternating least squares. 
In Section <ref> we develop several algorithms that rely on numerical 
tensor methods to solve high-dimensional PDEs. 
Specifically,  we discuss schemes based on implicit and explicit 
time integration, combined with dimensional splitting, 
alternating least squares and hierarchical singular value decomposition. 
In Section <ref> we demonstrate the effectiveness 
and computational efficiency of the proposed new algorithms 
in solving both an advection equation in six 
variables plus time and a linearized version of the Boltzmann equation.





§ TENSOR DECOMPOSITION OF HIGH-DIMENSIONAL FUNCTIONS



Consider a multivariate scalar field 
f( z): D⊆ℝ^N →ℝ. 
In this section we briefly review effective representations 
of f based on tensor methods. 
In particular, we discuss the canonical tensor decomposition 
and hierarchical tensor methods, e.g., tensor train and 
hierarchical Tucker expansions.



 §.§ Canonical Tensor Decomposition



The canonical tensor decomposition of the multivariate function
f( z): D⊆ℝ^N →ℝ is a 
series expansion of the form 

    f( z)≃∑_l=1^r ∏_k=1^NG^l_j(z_j),

where G^l_i(z_i): ℝ→ℝ are one-dimensional 
functions usually represented relative to a 
known basis {ϕ_1,...,ϕ_Q}, i.e., 

    G^l_i(z_i)=∑_k=1^Qβ_ik^l ϕ_k(z_i).

The quantity r in (<ref>) is 
called separation rank.

Although general/computable theorems relating a prescribed accuracy of the representation of f( z) to the value of the separation rank r are still lacking, there are cases where the 
expansion (<ref>) is exponentially 
more efficient than one would expect 
a priori <cit.>. 









  
Alternating Least Squares (ALS) Formulation
Development of robust and efficient algorithms to 
compute (<ref>) to any desired accuracy is 
still a relatively open question (see
<cit.> 
for recent progresses). Computing the tensor 
components G_k^l(z_k) usually relies on (greedy) 
optimization techniques, such as alternating least 
squares (ALS) <cit.> 
or regularized Newton methods <cit.>, which are only locally 
convergent <cit.>, i.e., the final result may depend 
on the initial condition of the algorithm. 

In the least-squares setting, the tensor 
components G^l_j(z_j) are computed by minimizing a norm of the residual,


    R( z)=f( z)-∑_l=1^r ∏_k=1^N G^l_j(z_j),
  
with respect to the rNQ degrees of freedom β_hj^n
i.e., 

    min_β_hj^nR( z).


Assuming f( z) to be periodic in the hyper-cube D=[-b,b]^N (b>0),
 we define the norm · as a standard L_2 norm

    R( z)_L_2^2 =   ∫_-b^b⋯∫_-b^b 
    R( z)^2 d z,
    
    =   ∫_-b^b⋯∫_-b^b 
    [f( z) -∑_n=1^r∏_k=1^N (∑_s=1^Q
    β_ks^nϕ_s(z_k))]^2   d z.


In the alternating least-squares (ALS) paradigm, 
we compute the minimizer of the residual (<ref>)
by splitting the non-convex optimization problem (<ref>) into 
a sequence of convex low-dimensional convex optimization
problems. To illustrate the method, let us define vectors 

    β_1 = (β^1_11, ..., β^1_1Q,...,β^r_11, ..., β^r_1Q)^T,   ⋯,  β_N=(β^1_N1, ..., β^1_NQ,...,β^r_N1, ..., β^r_NQ)^T.

Each β_i collects the degrees of freedom 
of all functions {G^1_i(z_i), ..., G^r_i(z_i)} depending 
on z_i. Next, the optimization 
problem (<ref>) is split into a sequence of 
convex optimization problems,

    min_β_1 R_L_2^2, 
    
      ⋯,   min_β_N R_L_2^2.

This sequence is not equivalent to the full problem (<ref>) because, 
in general, it does not allow one to compute the global minimizer of (<ref>) <cit.>. 


The Euler-Lagrange equations associated 
with (<ref>) are of the form[
Recall that minimizing the residual (<ref>)
with respect to β_hj^n 
is equivalent to impose orthogonality 
relative to the space spanned by the functions

    ϕ_h(z_j)∏^N_k=1
    k≠ jG^n_k(z_k).

]

    A_j β_j = g_j     j=1,…,N,


where

    A_j=[
    [ A^11_j11⋯ A^11_j1Q A^12_j11⋯ A^12_j1Q                  ⋯ A^1r_j11⋯ A^1r_j1Q;                  ⋮                  ⋮                  ⋮                  ⋮; A^11_jQ1⋯ A^11_jQQ A^12_jQ1⋯ A^12_jQQ                  ⋯ A^1r_jQ1⋯ A^1r_jQQ;                  ⋮                  ⋮                  ⋮                  ⋮; A^r1_j11⋯ A^r1_j1Q A^r2_j11⋯ A^r2_j1Q                  ⋯ A^rr_j11⋯ A^rr_j1Q;                  ⋮                  ⋮                  ⋮                  ⋮; A^r1_jQ1⋯ A^r1_jQQ A^r2_jQ1⋯ A^r2_jQQ                  ⋯ A^rr_jQ1⋯ A^rr_jQQ ]],     g_j=[
    [ g_j1^1;      ⋮; g_jQ^1;      ⋮; g_j1^r;      ⋮; g_jQ^r;        ]],

and 

    A^ln_jhs=∫_-b^b ϕ_h(z_j) ϕ_s(z_j) dz_j 
    ∏_k=1
    k≠ j^N
    ∫_-b^b G_k^l(z_k)G_k^n(z_k) dz_k,


    g_jh^n = ∫_-b^b⋯∫_-b^b f( x) ϕ_h(x_j)∏_k=1
    k≠ j^N
    G_k^n(x_k) dz_1⋯dz_N.


The matrices A_j are symmetric, positive 
definite and of size rQ× rQ. 




  
Convergence of the ALS Algorithm 
The ALS algorithm described above is  
an alternating optimization scheme, i.e.,  
a nonlinear block Gauss–Seidel 
method (<cit.>, 7.4). 
There is a well–developed local 
convergence theory for this type of methods <cit.>. 
In particular, it can be shown that ALS is locally 
equivalent to the linear block Gauss–Seidel 
iteration applied to the Hessian matrix. 
This implies that ALS is linearly convergent in the iteration 
number <cit.>, 
provided that the Hessian of the residual 
is positive definite (except on a trivial null 
space associated with the scaling non-uniqueness 
of the canonical tensor decomposition). The last assumption 
may not be always satisfied. Therefore, convergence 
of the ALS algorithm cannot be granted in general.
Another potential issue of the ALS algorithm is 
the poor conditioning of the matrices A_j 
in (<ref>), which can addressed by 
regularization <cit.>.

The canonical tensor decomposition 
(<ref>) in N dimensions has relatively 
small memory requirements. In fact, the number of degrees of freedom 
that we need to store is rN  Q, where r is 
the separation  rank, and Q is the number of degrees of 
freedom employed in each tensor component G_k^l(z_k).
Despite the relatively low-memory requirements, it is often desirable 
to employ scalable parallel versions the ALS algorithm 
<cit.> to compute the canonical 
tensor expansion (<ref>) 
(see Appendix <ref>)





 §.§ Hierarchical Tensor Methods


Hierarchical Tensor methods <cit.>  
were introduced to mitigate the dimensionality problem in the  
core tensor of the classical Tucker decomposition <cit.>. 
A key idea is to perform a sequence of 
Schmidt decompositions (or multivariate SVDs <cit.>) until the 
approximation problem is reduced to a product 
of one-dimensional functions/vectors. 
To illustrate the method in a simple way, 
consider a six-dimensional function
f( z)= f(z_1,…,z_6). 

We first split the variables as (z_1,z_2,z_3) and 
(z_4,z_5,z_6) through one Schmidt 
decomposition <cit.> as

    f( z)=    ∑_i_7,i_8=1^r A^{1}_i_7i_8 T^{1,2,3}_i_7
    (z_1,z_2,z_3)T^{4,5,6}_i_8(z_4,z_5,z_6).

Then we decompose T^{1,2,3}_i_7(z_1,z_2,z_3) and 
T^{4,5,6}_i_8(z_4,z_5,z_6) further by 
additional Schmidt expansions to obtain
 
    f( z)=   ∑_i_7,i_8=1^r A^{1}_i_7i_8∑_i_1,i_9=1^r  A^{2}_i_7i_1i_9
     T^{1}_i_1(z_1)
     T^{2,3}_i_9(z_2,z_3)
    ∑_i_4,i_10=1^r A^{3}_i_8 i_4i_10
    T^{4}_i_4(z_4)T^{5,6}_i_10(z_5,z_6),
       =∑_i_1,⋯, i_6=1^r
    C_i_1⋯ i_6
    T^{1}_i_1(z_1)T^{2}_i_2(z_2)T^{3}_i_3(z_3)
    T^{4}_i_4(z_4)T^{5}_i_5(z_5)T^{6}_i_6(z_6).

The 6-dimensional core tensor[A diagonalization 
of the core tensor in (<ref>) would minimize the number of 
terms in the series expansion. Unfortunately, this is impossible for a tensor with 
dimension larger than 2 (see, e.g., <cit.>)
and, for complex tensors, <cit.>). A closer look at the canonical tensor 
decomposition (<ref>) reveals that such an expansion is 
in the form of a fully diagonal high-order Schmidt decomposition, i.e., 

    f(x_1,…,x_N)=∑_l=1^r C_l⋯ l
    G^l_1(x_1)⋯ G^l_N(x_N).

The fact that diagonalization of C_l_1⋯  l_6 is 
impossible in dimension larger than 2 implies that it is impossible 
to compute the canonical tensor decomposition of f by standard 
linear algebra techniques.] 

is explicitly obtained as


    C_i_1⋯ i_6 =
    ∑_i_7,i_8=1^r A^{1}_i_7i_8∑_i_9=1^r  A^{2}_i_7i_1i_9 A^{4}_i_9i_2i_3∑_i_10=1^r A^{3}_i_8 i_4i_10A^{5}_i_10i_5i_6.

The procedure just described forms the foundation of 
hierarchical tensor methods. The key element is that 
the core tensor C_i_1⋯ i_6 resulting from 
this procedure is always factored as 
a product of at most three-dimensional matrices. 
This is also true in arbitrary dimensions. 
The tensor components T^{k}_i_j and the factors of the  
core tensor (<ref>) can be effectively 
computed by employing hierarchical singular value decomposition  <cit.>.
Alternatively, one can use an optimization framework 
that leverages recursive subspace factorizations 
<cit.>. 
If we single out one variable at the time and perform a 
sequential Schmidt decomposition of the remaining 
variables we obtain the so-called tensor-train (TT)
decomposition <cit.>. 
Both tensor train and hierarchical tensor expansions
can be conveniently visualized by  
graphs (see Figure <ref>). This is done 
by adopting the following standard rules: i)  
a node in a graph represents a tensor in as many variables as 
the number of the edges connected to it, ii) 
connecting two tensors by an edge represents 
a tensor contraction over the index 
associated with a certain variable. 



Efficient algorithms to perform basic operations 
between hierarchical tensors, such as addition, 
orthogonalization, rank reduction, scalar products, 
multiplication, and linear transformations, are discussed in 
<cit.>.
Parallel implementations of such algorithms were recently proposed in 
<cit.>.

Methods for reducing the computational cost of tensor 
trains are discussed in <cit.>. 
Applications to the Vlasov kinetic equation can be found in 
<cit.>.










§ NUMERICAL APPROXIMATION OF  HIGH-DIMENSIONAL PDES



In this section we develop efficient numerical methods to integrate 
high-dimensional linear PDEs of the form 

    ∂ f( z,t)/∂ t = L( z) f( z,t),

where L( z) is a linear operator. 
The algorithms we propose are based on numerical tensor 
methods and appropriate rank-reduction techniques, such as 
hierarchical singular value decomposition <cit.> and 
alternating least squares <cit.>.

To begin with, we assume that L( z) is a separable 
linear operator with separation rank r_L, i.e.,  an operator of the form


    L( z) = ∑_l=1^r_Lα_l L_1^l(z_1)⋯ L_N^l(z_N).


For each i and j, L_i^j(z_i) is a linear operator acting only on variable z_i.

As an example,  the operator

    L(z_1,z_2)= z^2_1∂/∂ z_1-sin(z_2)∂^2 /∂ z_1∂ z_2.

is separable with separation rank r_L=2 in dimension N=2. It can be written in the form (<ref>) if we set  
α_1=α_2=1 and 

    L_1^1(z_1)=z_1^2∂/∂ z_1,    L_2^1(z_2)=1,   
    L_1^2(z_1)=∂/∂ z_1,    L_2^2(z_2)=-sin(z_2)∂/∂ z_2.







 §.§ Tensor Methods with Implicit Time Stepping


Let us discretize the PDE in (<ref>) in time by using the 
Crank-Nicolson method. To this end, consider an evenly-spaced 
grid t_n= nΔ t (n=0,1,...) with time step Δ t.  This yields 

    [I-Δ t/2 L( z) ]f_n+1( z)=
    [I+Δ t/2 L( z) ]f_n ( z)+ Δ tτ_n+1( z),

where I is the identity matrix, f_n( z)=f( z,t_n), and τ_n+1( z) is the 
local truncation error of the Crank-Nicolson method 
at time t_n+1 <cit.>.
Rewriting this in a more compact notation yields <cit.>

    A( z) f_n+1( z)=B( z) f_n( z)+Δ t τ_n+1( z),

where  

    A = I-Δ t/2 L( z),   and   
    B = I+Δ t/2 L( z).


It follows from the definition of L( z) in (<ref>) that both A( z) and B( z) are separable operators 
of the form  

    A= ∑_q=0^r_Lη_q E^q_1(z_1)⋯ E^q_N(z_N),      
    B= ∑_q=0^r_Lζ_q E^q_1(z_1)⋯ E^q_N(a_N),

where η_0=ζ_0=1, E^0_j(z_j)=1 (j=1,...,N), 

    η_j=-Δ t/2,    ζ_j=-η_j,     E^q_j(z_j)=L^q_j(z_j),      q=1,...,r_L,    j=1,...,N.


A substitution of the canonical tensor  
decomposition[Recall that the 
functions G_k^l(z_k,t_n) 
are in the form 

    G_k^l(z_k,t_n)=∑_s=1^Q 
    β^l_ks(t_n)ϕ_s(z_k),

where β^l_ks(t_n) (l=1,...,r, k=1,...,N, s=1,...,Q) 
are the degrees of freedom. ]


    f̂_n+1( z) = ∑_l=1^r ∏_k=1^N G^l_k(z_k,t_n+1)

into (<ref>) yields the residual 

    R( z,t_n+1) = A( z)f̂_n+1( z)-B( z)f̂_n( z),

in which the local truncation error Δ t τ_n+1 is embedded into R( z,t_n+1).
Next, we minimize the residual using the alternating least squares
algorithm, described in Section <ref>, to obtain the solution at time t_n+1.

In particular, we look for a minimizer of (<ref>) 
computed in a parsimonious way. 

The key idea again is to split the optimization problem


    min_β_ks^l(t_n+1)R( z,t_n+1)_L_2^2


into a sequence of optimization problems of smaller
dimension, which are solved sequentially and in parallel 
<cit.> (see Appendix <ref>). 
To this end, we define

    β_k(t_n)=[β^1_k1(t_n), ..., β^1_kQ(t_n),...,
    β^r_k1(t_n), ..., β^r_kQ(t_n)]^T      k=1,...,N.

The vector β_k(t_n) collects 
the degrees of freedom representing the solution 
functional along z_k at time t_n, i.e., 
the set of functions {G_k^1(z_k,t_n),...,G^r_k(z_k,t_n)}.

Minimization of (<ref>) with respect to 
 independent variations of β_k(t_n+1)
yields a sequence of convex optimization problems (see Figure <ref>)

    min_β_1(t_n+1)R( z,t_n+1)_L_2^2,   min_β_2(t_n+1)R( z,t_n+1)_L_2^2,   ⋯,  min_β_N(t_n+1)R( z,t_n+1)_L_2^2.



This set of equations defines the alternating least-squares (ALS) method. 

The Euler-Lagrange equations, which identify stationary points 
of (<ref>), are linear systems in the form

    M^L_q β_q(t_n+1) =  M^R_q β_q(t_n)     q=1,...,N

where

    M^L_q = ∑_e,z=0^r_L K^L_ez[
    
    
    
     N
    ◯
    k=1
    k≠ q
    β̂_k(t_n+1)^T  
    E^ez_k β̂_k(t_n+1)
    ]^T⊗[ E^ez_q]^T,


    M^R_q = ∑_e,z=0^r_LK^R_ez[ 
    
    
    
     N
    ◯
    k=1
    k≠ q
    β̂_k(t_n)^T 
     E^ez_k β̂_k(t_n+1)
    ]^T⊗[ E^ez_q]^T,

and

    [ E^ez_q]_sh=∫_-b^b E_q^e(a)ϕ_s(a) E_q^z(a)ϕ_h(a) da.

In (<ref>)–(<ref>), 
◯ denotes the Hadamard matrix 
product, ⊗ is the Kronecker matrix product, 
β̂_k is the matricization of β_k, i.e., 

    β̂_k(t_n)=
    [[  β^1_k1(t_n)            ⋯  β^r_k1(t_n);            ⋮            ⋱            ⋮;  β^1_kQ(t_n)            ⋯ β^r_kQ (t_n) ]],

E^ez_q are Q× Q matrices, 
and  K^L_ez and K^R_ez are entries of the 
matrices

    K^L = ηη^T,  and   K^R = ηζ^T,

where η and ζ are column vectors with entries η_q and ζ_q defined in (<ref>).
The ALS algorithm, combined with 
canonical tensor representations, effectively reduces evaluation of N-dimensional integrals to evaluation of the sum of products of one-dimensional integrals.



  
Summary of the Algorithm
Often many of the integrals in (<ref>) only differ by a prefactor and in that
case the number of unique integrals can be reduced to a very small number
(4 in the case of the Boltzmann-BGK
equation discussed in Section <ref>). 
Thus, it is convenient to precompute a map between such 
integrals and any entry in (<ref>).
Such a map allows us to rapidly compute each 
matrix E_q^ez  
in (<ref>) and (<ref>), and therefore to efficiently build the matrix system. 

Next, we compute the canonical tensor 
decomposition of the initial condition f( z,0)
by applying the methods described in Section <ref>.
This yields the set of vectors 
{β_1(t_0),...,β_N(t_0)}, from which the matrices 
M^L_1 and M^R_1 in (<ref>) and 
(<ref>) are built. To this end, we need an initial guess 
for {β_1(t_1),...,β_N(t_1)} given, e.g., by {β_1(t_0),...,β_N(t_0)} or its small random perturbation. 
With M^L_1, M^R_1 in place, 
we solve the linear system in (<ref>) and 
update β_1(t_1). Finally, we recompute  M^L_2 and M^R_2  
(with the updated β_1(t_1)) and solve for β_2(t_1). 
This process is repeated for q=3,…,N and iterated among 
all the variables until convergence (see Figure <ref>). 
The ALS iterations are said to converge if 

    β^(j)_k(t_n)-β^(j-1)_k(t_n)_2^2≤ϵ_Tol,      j=1,…,N_max,

where j here denotes the iteration number whose maximum value is N_max, and ϵ_Tol is a prescribed tolerance on the increment of such iterations. 

Parallel versions of the ALS algorithm associated with canonical 
tensor decompositions were recently proposed in <cit.> 
(see also Appendix <ref>).

We emphasize that the use of a hierarchical tensor series 
instead of (<ref>) (see Section (<ref>)) also allows one to compute the unknown degrees of the expansion by using parallel ALS algorithms. This question was recently addressed in <cit.>.




 §.§ Tensor Methods with Explicit Time Stepping


Let us discretize PDE (<ref>) in time by using any explicit 
time stepping scheme, for example the second-order Adams-Bashforth 
scheme <cit.>


    f_n+2( z)=f_n+1( z)+Δ t/2 L( z)(3f_n+1( z)-f_n( z)) +Δ t τ_n+2( z).


In this setting, the only operations needed to 
compute f_n+2 with tensor methods are: i) addition, 
ii) application of a (separable) linear operator L, and iii) rank 
reduction[On the other hand, the use of
implicit time-discretization schemes requires one to  
 solve linear systems with tensor methods <cit.>.], 
the last operation being the most important among the three
<cit.>. 
Rank reduction can be achieved, e.g., by hierarchical singular value 
decomposition <cit.> or by 
alternating least squares <cit.>.

From a computational perspective, it is useful to 
split the sequence of tensor operations yielding f_n+2 in (<ref>) 
into simple tensor operations followed by rank reduction. For example, 
 L(3f_n+1-f_n) might be evaluated as follows: i) compute 
w_n+1=f_n+1-f_n; ii) perform rank reduction on w_n+1; iii) compute 
L w_n+1 and perform rank reduction again.  
This allows one to minimize the memory requirements and overall 
computational cost.[Recall that adding two tensors 
with rank r_1 and r_2 usually yields a tensor of rank r_1+r_2.]
Unfortunately, splitting tensor operations into sequences of tensor operations 
followed by rank reduction can produce severe cancellation 
errors. In some cases this problem can 
be mitigated. For example, an efficient and robust algorithm, which 
allows us to split sums and rank reduction operations, was recently proposed 
in <cit.> in the context of 
hierarchical Tucker formats. The algorithm leverages  
the block diagonal structure arising from addition of
hierarchical Tucker tensors. 




§ NUMERICAL RESULTS


In this section, we demonstrate the accuracy and 
effectiveness of the numerical tensor methods discussed above. Two case studies are considered: 
an advection equation in six dimensions plus time and the Boltzmann-BGK equation.




 §.§ Advection Equation

Consider an initial-value problem

    ∂ f/∂ t+∑_j=1^N(∑_k=1^N 
    C_jkz_k)∂ f/∂ z_j=0,
         f(z_1,…,z_N,0) = f_0(z_1,…,z_N),


where N is the number of independent variables forming the coordinate vector z = (z_1, …, z_N)^⊤, and C_jk are components of a given matrix of coefficients C. 
The analytical solution to this initial-value problem is computed with the method of 
characteristics as <cit.>

    f( z,t)=f_0 ( e^-t C  z).

We choose the initial condition f_0 to be a fully separable product of exponentials,

    f_0( z)=e^- z ^2_2.


The time dynamics of the solution (<ref>)-(<ref>)
is entirely determined by the matrix C.
In particular, if C has eigenvalues with 
positive real part then the matrix exponential 
exp(-t C) is a contraction map. 
In this case, f( z,t)→ 1 as t→∞,
at each point z∈ℝ^N. 

Figure <ref> exhibits the analytical solution (<ref>) 
generated by a 2× 2 matrix C, whose eigenvalues are complex conjugate with positive real part.



The side length of the hypercube cube that encloses 
any level set of the solution at time t depends on 
the number of dimensions N. In particular, the side-length 
of the hypercube that encloses 
the set { z∈ℝ^N |  f( z,t)≥ϵ}
is given by 

    1/2√(-log(ϵ)/ N λ_min)≤ b 
    ≤1/2√(-log(ϵ)/λ_min).

where λ_min is the smallest eigenvalue of the matrix 
exp( t  C)^⊤exp(t  C).
Equation (<ref>) can be written in the form 
(<ref>), with separable operator L (separation rank N^2). 
Specifically, the one-dimensional operators L^l_j(z_j) appearing 
in (<ref>) are explicitly defined in Table <ref>.



We used tensor methods, both canonical tensor decomposition and 
hierarchical tensor methods, and explicit time stepping (Section <ref>) to solve numerically (<ref>) in the 
periodic hypercube[Such domain is chosen large 
enough to accommodate periodic (zero) boundary 
conditions in the integration period of interest.] [-60,60]^N.  Each tensor component was discretized by using a Fourier spectral expansion with 600 nodes in each variable (e.g., basis functions in (<ref>) or (<ref>)). 

The accuracy of the numerical solution was quantified in terms of the time-dependent relative error

    ϵ_m(t) = |f( z^*,t)-f̂( z^*,t)/f( z^*,t)|,
  
where f is the analytical solution (<ref>), f̂ 
is the numerical solution obtained by using the 
canonical or hierarchical tensor methods with separation rank r. 



Figure <ref> shows the relative pointwise error (<ref>), computed at z^*=(h,h,…,h) with h= 0.698835274542439, for different separation ranks r and  
different number of independent variables N. 
As expected, the accuracy of the numerical solution increases with the separation rank r. 
Also, the relative error increases with the number of dimensions N.

It is worthwhile emphasizing that the rank reduction in 
canonical tensor decomposition at each time step is based 
on a randomized algorithm that requires 
initialization at each time-step. This means that 
results of simulations with the same nominal 
parameters may be different. On the other hand, tensor
methods based on a hierarchical singular 
value decomposition, such as the hierarchical Tucker 
decomposition <cit.>, do 
not suffer from this issue.
The separation rank of both canonical and hierarchical tensor 
methods is computed adaptively up 
to the maximum value r_max specified 
in the legend of Figure <ref>. 

In two dimensions, the CP-ALS algorithm 
yields the similar error plots for 
r_max=8 and r_max=12. This is because r < 8 in both cases, and 
throughout the simulation up to t=1. 
The difference between these results is due to the 
random initialization required by the ALS algorithm 
at each time step. 





 §.§ Boltzmann-BGK Equation


In this section we develop an accurate ALS-based algorithm  to solve  a linearized BGK equation (see Appendix <ref> for details of its derivation),


    ∂ f/∂ t+ v ·∇_ x f = 
    ν I,      I = f_eq( v)-f.


Here, f( z,t) is a probability density function is six phase variables plus 
time; the coordinate vector z = ( x,  v)^⊤ is composed of three spatial dimensions x = (x_1,x_2,x_3)^⊤ and three components of the velocity vector v = (v_1,v_2,v_3)^⊤; f_eq denotes a locally Maxwellian distribution,

    f_eq( v) = ρ/(2π R T)^3/2exp(- v_2^2/2R T),


with uniform gas density ρ, temperature T, and velocity v; R is the gas constant; and ν = κρ T^1-μ with positive constants κ and μ. The solution to the linearized Boltzmann-BGK equation (<ref>) is also
computable analytically, which provides us with a benchmark solution 
to check the accuracy of our algorithms. 

As before we start by rewriting this equation 
in an operator form, 

    ∂ f/∂ t = L f + C

where 

    L( x, v) = - v·∇_ x- ν I,    
    C( v) = ν f_eq( v).

Note that L is a separable linear 
operator with separation rank r_L=4, which can be rewritten in the form 

    L( z)= ∑_q=1^4α_q L^q_1(z_1)⋯ L_6^q(z_6),

for suitable one-dimensional linear operators L^q_j(a_j) 
defined in Table <ref>. Similarly, C( v) in (<ref>) is separated as 

    C( v) = C_1(v_1)C_2(v_2)C_3(v_3),   where   
    C_i(v_i)=(νρ)^1/3/(2π R T)^1/2exp(-v_i^2/2R T).





From a numerical viewpoint, the solution to (<ref>) can be 
represented by using any of the tensor series expansion
we discussed in Section <ref>. In particular, 
hereafter we develop an algorithm based on canonical tensor tensor 
decomposition, alternating least squares, and implicit time 
stepping (see Section (<ref>)).
To this end, we discretize (<ref>) in time with 
the Crank-Nicolson method to obtain

    [I-Δ t/2 L ]f_n+1=
    [I+Δ t/2 L ]f_n+Δ t C + Δ tτ_n,

where f_n( z)= f( z,t_n) and τ_n it the local truncation 
error of the Crank-Nicholson method at t=t_n (<cit.>, p. 499). 
This equation can be written in a compact notation as 

    Af_n+1=Bf_n + Δ t C +Δ t τ_n,

where  A and B are separable operators in the form (<ref>), 
where all quantities are defined in Table <ref>.



A substitution of the canonical tensor  
decomposition[Recall that the 
functions G_k^l(z_k,t_n) 
are in the form 

    G_k^l(z_k,t_n)=∑_s=1^Q 
    β^l_ks(t_n)ϕ_s(z_k),

β^l_ks(t_n) (l=1,...,r, k=1,...,6, s=1,...,Q) 
being the degrees of freedom. ]


    f̂_n+1 = ∑_l=1^r ∏_k=1^6 G^l_k(z_k,t_n+1)

into equation (<ref>) yields the residual 

    R( z,t_n+1) = A( z)f̂_n+1( z)-B( z)f̂_n( z) + Δ t C(z_4,z_5,z_6),

which can be minimized by using the alternating least squares method, as
we described in Section <ref> 
to obtain the solution at time t_n+1.






  
Nonlinear Boltzmann-BGK Model

The numerical tensor methods we discussed in Section <ref> and Section <ref> can be extended 
to compute the numerical solution of the fully nonlinear Boltzmann-BGK model

    ∂ f/∂ t + v·∇_ x f = 
     ν( x,t)(f_eq( x, v,t)-f( x, v,t)).

Here, f_eq( x, v,t) is the equilibrium distribution (<ref>), 
while the collision frequency ν( x,t) can can be expressed, e.g., 
as in  (<ref>). From a mathematical viewpoint both quantities 
f_eq and ν are nonlinear functionals of the PDF f( x, v,t) 
(see equations (<ref>)-(<ref>)). Therefore (<ref>) is  
an advection equation driven driven by a nonlinear functional 
of f( x, v,t). The evaluation of the integrals appearing 
in (<ref>)-(<ref>) poses no great challenge if we represent f as  
a canonical tensor (<ref>) or as a hierarchical Tucker tensor
(<ref>). Indeed, such integrals can be factored as sums of 
product of one-dimensional integrals in a tensor representation. Moreover, 
computing the inverse of ρ( x,t) in 
(<ref>) and (<ref>) is relatively simple in a collocation 
setting. In fact, we just need to evaluate the full tensor ρ( x,t), on 
a three-dimensional tensor product grid and 
compute  1/ρ( x,t) pointwise. If needed, we could then compute 
the tensor decomposition of 1/ρ(x,t), e.g., by using the canonical tensor representations and the ALS algorithm we discussed in Section <ref>.
This allows us to evaluate the BGK collision operator  and represent 
it in a tensor series expansion. The solution to the Boltzmann BGK 
equation can be then computed by operator splitting methods <cit.>.







  
Initial and Boundary conditions
To validate the proposed alternating least squares algorithm we 
set the initial condition to be either coincident with the homogeneous
Maxwell-Boltzmann distribution (<ref>), i.e.,


    f ( x, v,0) =
      ρ/(2 π R T)^3/2exp(- v_2^2/2 R T)

or a slight perturbation of it. 
The PDF (<ref>) is obviously separable with separation rank one. 
Moreover, we set the computational domain to be a
the periodic hyperrectangle [-b_x,b_x]^3×[-b_v,b_v]^3.
It is essential that such domain is chosen large
enough to accommodate the support of the PDF at all times, thus
preventing physically unrealistic correlations between the distribution function. 
More generally, it is possible to enforce other types of boundary conditions 
by using appropriate trial functions, or by a mixed approach
<cit.> in the case of 
Maxwell boundaries.

Unless otherwise stated, all simulation parameters  
are set as in Table <ref>. Such setting corresponds 
to the problem of computing the dynamics of Argon within 
the periodic hyperrectangle [-b_x,b_x]^3× [-b_v,b_v]^3. 
We are interested in testing our algorithms in two different regimes: 
i) steady state, and ii) relaxation to statistical equilibrium. 






  §.§.§ Canonical Tensor Decomposition of the Maxwell-Boltzmann Distribution

We begin with a preliminary convergence study of the canonical tensor 
decomposition (<ref>) applied to the 
one-particle Maxwell distribution f_eq( v). Such study 
allows us to identify the size of the hyperrectangle that 
includes the support of the Maxwellian equilibrium  distribution. 
This is done in Figure <ref>, where we plot 
the tensor expansion of f_eq( v) as function 
of the non-dimensional molecular speed v_2/√(R_S T), 
for different number of basis functions Q in (<ref>), and for
various values of the dimensionless velocity domain size b_v/√(R T). 
It is seen that for small values of b_v/√(R T), a small number of 
basis functions Q is sufficient to capture the exact equilibrium 
distribution. At the same however, small values of
b_v/√(R T) result in distribution functions with 
truncated tails, as can be seen in Figures <ref> (a), <ref>
(b), and <ref> (e). 

Figures <ref> (c) and <ref> (g)
demonstrate that the dimensionless velocity domain size b_v/√(R  T) = 5 
is sufficient to avoid truncation of the tail of the distribution, 
and that a series expansion with Q=11 in each variable is sufficient 
to accurately capture the equilibrium distribution. This justifies 
the choice of parameters in Table (<ref>).



A more detailed analysis of the effects of the truncation order Q
and the dimensionless velocity domain size b_v/√(R T) 
on the approximation error in f_eq is done in Figure <ref>(a), 
where we plot the Normalized Mean Absolute Error (NMAE) 
versus b_v/√(R T) and Q. The NMAE between two 
vectors X and Y of 
size N is defined as: 

    NMAE 
    = 1/N X - Y_1/max( X)-min( Y),

where ·_1 is the standard vector 1-norm. 
In Figure  <ref>(b) we plot the NMAE as a function of 
Q √(R T)/b_v. This results in a collapse of the data which suggests 
that if one wants to double the non-dimensional velocity domain size while
maintaining the same accuracy, the series truncation order Q needs 
to be doubled as well. In addition, it can be seen that 
starting at Q √(R T)/b_v≈ 1 the error decays 
rapidly with increasing Q √(R T)/b_v.





  §.§.§ Steady State Simulations

We first test the the Boltzmann-BGK solver we developed 
in Section (<ref>) on an initial value problem 
where the initial condition is set to be the Maxwell-Boltzmann
distribution f_eq( v) (see equation (<ref>)). 
A properly working algorithm should keep such equilibrium distribution 
unchanged as the simulation progresses.
In Figure <ref> we plot the 
distribution of molecular velocities f( v,t) as function of the molecular 
speed v_2ξ⃗ at various dimensionless 
times t/τ_R, where τ_R is the relaxation time in Table <ref>. 
It is seen that the alternating least squares algorithm we developed in 
Section (<ref>) has a stable 
fixed point at the Maxwell-Boltzmann equilibrium 
distribution f_eq. This is an important test, as 
convergence of alternating least square is, in general, not granted 
for arbitrary residuals (see, e.g., <cit.>).



We also computed the moments of the distribution 
function f( x, v,t) at different times to verify 
whether ALS iterations preserve  (density), 
momentum (velocity), and kinetic energy (temperature). 
In particular, we computed  

    < ρ>(t)
       = 
      1/b_x^3∫_-b_x^b_x⋯∫_-b_v^b_v
      f( x, v,t)d x d  v 
    <  u >(t)
       =
      1/< ρ> b_r^3∫_-b_x^b_x⋯∫_-b_v^b_v v
      f ( x, v,t)d x d  v,    and 
    < T>(t)
       =
      1/< ρ>  b_x^3 R_S∫_-b_x^b_x⋯∫_-b_v^b_v v - <  u >_2^2
      f ( x, v,t)d x d  v

In Figure <ref> we plot(<ref>)-(<ref>) versus 
time. It is seen that such quantities are indeed constants, i.e., 
ALS iterations preserve the average density, velocity and temperature.








  §.§.§ Relaxation to Statistical Equilibrium

In this section we consider relaxation statistical equilibrium in the 
linearized Boltzmann-BGK model (<ref>). This  
allows us to study the accuracy and computational efficiency of the 
proposed ALS algorithm in transient dynamics. To this end, 
we consider the initial condition 

    f( x, v,0)
    = f_eq(1 + ϵ∏_k=1^3cos(2 πx_k/b_x))

with ϵ = 0.3. Such initial condition is a slight 
perturbation of the Maxwell-Boltzmann equilibrium distribution 
(see Figure <ref>(a)). The simulation results we obtain 
are shown in Figures <ref>(b)-(d), where we plot 
the time evolution of the marginal PDF f( v,t) versus the 
dimensionless time t/τ_R. It is seen that 
the initial condition (<ref>) is in the basin of
attraction of f_eq, i.e., the ALS algorithm sends 
f( x, v,0) into f_eq( v) after 
a transient. Moreover, f( v,t) is attracted by f_eq( v) 
exponentially fast in time at a rate τ_R = 1/ν 
(see Figure <ref>), consistently with 
results of perturbation analysis.


In Figure <ref> we demonstrate that the parallel ALS algorithm we developed conserves the average density of particle, velocity (momentum) and temperature throughout the transient that yields statistical equilibrium.   


In Figure <ref> we plot the the Normalized Mean Absolute Error
(NMAE) between the Maxwell-Boltzmann distribution f_eq( v) 
and the distribution we computed numerically by using the proposed 
ALS algorithm. Specifically we plot NMAE versus the dimensionless 
time t/τ_R for different time-steps Δ t we employed
in our simulations.  

The goal of such study is to assess the sensitivity 
of parallel ALS iterations on the time step Δ t, and in turn on 
the total number of steps for a fixed integration time T=5τ_R. 
Results of Figure <ref> demonstrate that the decay of the 
distribution function f( x, v,t) to the its equilibrium state 
is indeed exponential. Moreover, we see that 
the ALS algorithm is robust to Δ t, which implies that 
the such algorithm is not sensitive to the total number of time steps 
we consider within a fixed integration time. 



Finally, we'd like to address scalability of the parallel ALS algorithm 
we developed in  (<ref>)., i.e., 
performance relative to the number of degrees of freedom 
and number of CPUs. 

To this end, we analyze here a few test cases we ran on a small workstation 
with 4 CPUs. The initial condition is set as f_eq in all cases 
and we integrated the Boltzmann-BGK model (<ref>) for 
1000 steps, with Δ t= 0.01 τ_R. The separation rank of the 
solution PDF is set to r=1 (see equation (<ref>)). 

In Figure <ref> (a) we plot the wall time t_wall
versus the total number of degrees of freedom of the ALS algorithm, 
i.e. 6Qr. 

The plot shows that a typical simulation of the Boltzmann-BGK model 
in 6 dimensions with r = 1 and Q=11 takes about
t_wall≈ 80min on a single core to complete 
a simulation up t = 10. 
With such value of Q there is an error of only 0.5% compared to the
analytical Maxwell-Boltzmann distribution function. 
The total wall time scales with the power 5/2 relative to 
number of degrees of freedom 6Qr.
In Figure <ref> we show that the total wall time scales we obtain with the proposed parallel ALS algorithm scales almost linearly 
with the number of CPUs (power -4/5). The leveling
off of performance for 4 cores can be attributed to 
the fact that the desktop computer we employed for our simulations 
is optimized for inter-core communication. We expect better
scaling performance of our ALS algorithm in computer clusters with 
infinity band inter-core communication.





§ SUMMARY


In this paper we presented new parallel algorithms 
to solve high-dimensional partial differential equations
based on numerical tensor methods.
In particular, we studied canonical and hierarchical tensor  
expansions combined with alternating least squares
and high-order singular value decomposition. Both implicit 
and explicit time integration methods were presented 
and discussed. We demonstrated the accuracy and 
computational efficiency of the proposed methods 
in simulating the transient dynamics generated by
a linear advection equation in six dimensions plus time, 
and the well-known Boltzmann-BGK equation.

We found that the algorithms we developed are extremely fast and
accurate, even in a naive Matlab implementation. Unlike direct
simulation Monte Carlo (DSMC), the high accuracy of the numerical 
tensor methods we propose makes it suitable for studying 
transient problems. 
On the other hand, given the nature of the numerical tensor discretization, 
implementation of the proposed algorithms to complex domains is not 
straightforward. 




§ ACKNOWLEDGEMENTS

We would like to thank N. Urien for useful discussions. D. Venturi 
was supported by the AFOSR grant FA9550-16-1-0092.






§ BOLTZMANN EQUATION AND ITS APPROXIMATIONS



In the classical kinetic theory of rarefied gas dynamics, 
the gas flow is described in terms of a non-negative density 
function f( x, v,t) which 
provides the number of gas particles at time t 
with velocity v∈ℝ^3 at position  
x ∈ℝ^3. 
The density function f satisfies the Boltzmann 
equation <cit.>.
In the absence of external forces, such 
equation can be  written as

    ∂ f/∂ t+ v ·∇_ x f = Q(f,f),
 
where Q(f,f) is the collision integral describing 
the effects of internal forces due to particle interactions. 
From a mathematical viewpoint the collision integral is a 
functional of the particle density. Its form depends on the 
microscopic dynamics. For example, for classical  
rarefied gas flows <cit.>,

    Q(f,f)( x, v,t)=∫_ℝ^3∫_𝕊^2
    B( v, v_1,ω) | f( x, v',t)f( x, v_1',t)-
    f( x, v,t)f( x, v_1,t) | dωd v_1.

In this expression, 

    v' = 1/2( v+ v_1+ v-  v_1 _2 w),     v_1' = 1/2( v+ v_1- v-  v_1 _2 w)

and { v, v_1} represent, respectively, 
the velocities of two particles before and after the 
collision, ω is a vector on the unit sphere 𝕊^2, 
and  B( v, v_1,ω) is the collision kernel. 
Such kernel is a non-negative function depending on the (Euclidean) 
norm v-v_1_2 and on the scattering angle θ 
between the relative velocities before and after the collision

    cos(θ) = ( v-  v_1)·ω/ v - v_1_2.

Thus, B( v, v_1,ω) becomes

    B( v, v_1,ω)= v -  v_1 _2
    σ( v- v_1_2,cos(θ)),

σ being the scattering cross section function <cit.>.
The collision operator (<ref>) satisfies a system of 
conservation laws in the form 

    ∫_ℝ^3 Q(f,f) ( x, v, t)ψ( v) d v=0,

where ψ( v) is either {1, v, v_2^2}.
This yields, respectively,  conservation of mass, momentum and energy. 
Moreover, Q(f,f) satisfies the Boltzmann H-theorem 

    ∫_ℝ^3 Q(f,f)( x, v, t)log(f( x, v,t)) 
    d v ≤ 0

Such theorem implies that any equilibrium distribution function, i.e.
any function f for which Q(f,f) = 0, has the form of a locally Maxwellian
distribution

    f_eq( x, v,t) = ρ( x,t)/(2π R T( x,t))^3/2exp(- u( x,t) -  v_2^2/2R T( x,t)),

where R is the gas constant, and ρ( x,t), u( x,t) and T( x,t) are the 
density, mean velocity and temperature of the gas, defined as 

    ρ( x,t) =    ∫_ℝ^3 f( x, v,t)d  v, 
     u ( x,t) =    1/ρ( x,t)∫_ℝ^3 v f( x, v,t) d  v,
     
    T( x,t) =    1/3ρ( x,t)∫_ℝ^3 u( x,t) -  v _2^2 f( x, v,t) d  v .

From a mathematical viewpoint, the Boltzmann equation (<ref>) is 
a nonlinear integro-differential equation for a positive scalar field in 
six dimensions plus time. By taking suitable averages over small volumes 
in position space, it can be shown that the Boltzmann equation 
is consistent with the compressible Euler's 
equations <cit.>, and the 
Navier-Stokes equations <cit.>.



  
BGK Approximation
The simplest collision operator satisfying the conservation law (<ref>)  
and the condition 

    Q(f,f)=0 ⇔ f( x, v,t)=f_eq( x, v,t)

(see equation (<ref>)) was proposed 
by Bhatnagar, Gross and Krook in <cit.>. The 
corresponding model, which is known as the BGK model, is defined 
by the linear relaxation operator

    Q(f,f)=ν( x,t) [f_eq( x, v,t)-f( x, v,t)]     ν( x,t)>0.

The field ν is usually set to be proportional to the density 
and the temperature of the gas <cit.>

    ν( x,t)=κρ( x,t)T( x,t)^1-μ.

It can be shown that the Boltzmann-BGK model (<ref>) 
converges to the correct Euler equations of incompressible fluid dynamics 
with the scaling x'=ϵ x, t'=ϵ t, and in the 
limit ϵ→ 0.  
However, the model does not converge to the correct 
Navier-Stokes limit. The main reason is that it predicts an 
unphysical Prandtl number <cit.>, which is 
larger that one obtained with the full collision 
operator (<ref>). 
The correct Navier-Stokes limit can be recovered by more 
sophisticated BGK-type models, e.g., the 
ellipsoidal statistical BGK model <cit.>. 
 
 
Following <cit.>, 
we introduce additional simplifications, namely we assume 
that  ν( x,t) =ν is constant and assume the equilibrium density, temperature and velocity 
to be homogeneous across the spatial domain. Assuming u( x,t)= 0, this yields an equilibrium distribution in (<ref>). The last assumption effectively decouples the BGK collision 
operator from the probability density function f( x, v,t). 
This, in turn, makes the BGK approximation linear, i.e., a six-dimensional PDE (<ref>).
Using both the isothermal assumption T( x,t)=T <cit.>, and
constant density assumption ρ( x,t)=ρ <cit.> means that our code can only be used to obtain
solutions which are small perturbations from the equilibrium distribution.
The term ν[f_eq( v)-f] in (<ref>) 
represents the relaxation time approximation of the Boltzmann 
equation, ν is the collision frequency, here assumed to be 
homogeneous.




§ PARALLELIZATION OF THE ALS ALGORITHM FOR THE BOLTZMANN-BGK EQUATION


In this appendix we describe the parallel alternating least squares 
algorithms we developed to solve the linearized Boltzmann-BGK equation
(<ref>).

Algorithm <ref> shows the main Alternating 
Least Square (ALS) routine while algorithm
<ref>, <ref>, <ref>, and
<ref> show the different subroutines. Each time step,
n, the value of β_Old is updated and then β_New is
determined in an iterative manner. Both the computations of the matrices M_k
and the vectors γ_k, as well as the calculation of the vectors
β_New,k can be performed independently of each other, and thus can
be easily parallelized. The current implementation assigns one core per dimension,
but further optimization could be possible by implementing a parallel method to
determine β_New <cit.>.



At each time step the CreateMatrixM subroutine
(algorithm <ref>) updates the different elements of matrix
M_k for every iteration until convergence has been reached. The algorithm
iterates over every element of the matrix and performs the multiplication and
summation of the elements of M̂. All the integrals were performed
analytically and the results are stored in the map Mk,l,z,s,q. 
Using this approach delivers significant savings in the 
computational cost of updating the matrix every iteration. 



The subroutine  (algorithm <ref>) updates
the vector γ_k. Similarly to the subroutine for the creation of matrix
M_k this algorithm iterates over all elements of the vector, and fills each
of them using the values of β_New and β_Old, and
pre-calculated maps. 








With the completion of the calculation of M_k and γ_k,
β_New can be updated using the  subroutine
(algorithm <ref>). Matrix M_k, depending on the rank
r_G, and series truncation Q, can become very large. In addition, it is
not known whether the matrix is invertible. This is why the system M_kβ_k = γ_k is solved using the iterative least square method, LSQR
<cit.>. For every iteration the value for
β_New,k which was determined in the previous iteration is used as
the initial guess for the algorithm.



Convergence of the system is calculated using the  subroutine
(algorithm <ref>). Because it is very computationally
expensive to calculate the full residual every iteration, instead the
convergence is determined using:

    ϵ_Conv 
    = 
      maxβ_New,k -
      β_Int,k/β_Int,k

and when ϵ_Conv < ϵ_Tol, the while loop in
algorithm <ref> is allowed to exit. In its currently implementation the
algorithm is able to reach convergence using a fixed Rank G, r_G = 2.
However, a future implementation could implement a dynamic Rank G for improved
accuracy. Every iteration the new value of
β_New is calculated according to:

    β_Newβ_Int + β_New - β_Int/Δβ

with Δβ = 4. To find the location where the residual reaches a
minimum, it is important to gradually approach this location and not update
β_New too aggressively by taking a smaller value of Δβ.
This causes overshooting of the minimum and the prevents convergence of
β_New and thus minimization of  the residual.






































































































































































































































































































































































































































































































































§ REFERENCES

100
url<#>1urlprefixURL href#1#2#2 #1#1

cercignani1988
C. Cercignani, The Boltzmann Equation and Its Applications, Springer New
  York, 1988.
http://dx.doi.org/10.1007/978-1-4612-1039-9
  doi:10.1007/978-1-4612-1039-9.

moss1995
F. Moss, P. V. E. McClintock (Eds.), Noise in nonlinear dynamical systems,
  Cambridge University Press, 1989.
http://dx.doi.org/10.1017/cbo9780511897832
  doi:10.1017/cbo9780511897832.

sobczyk2001
K. Sobczyk, Stochastic Differential Equations, Springer Netherlands, 1991.
http://dx.doi.org/10.1007/978-94-011-3712-6
  doi:10.1007/978-94-011-3712-6.

venturi2012
D. Venturi, T. P. Sapsis, H. Cho, G. E. Karniadakis, A computable evolution
  equation for the joint response-excitation probability density function of
  stochastic dynamical systems, Proceedings of the Royal Society A:
  Mathematical, Physical and Engineering Sciences 468 (2139) (2011) 759–783.
http://dx.doi.org/10.1098/rspa.2011.0186
  doi:10.1098/rspa.2011.0186.

li2009
J. Li, J. Chen, Stochastic Dynamics of Structures, John Wiley & Sons, Ltd,
  2009.
http://dx.doi.org/10.1002/9780470824269
  doi:10.1002/9780470824269.

shlesinger1998
M. F. Shlesinger, T. Swean, Stochastically Excited Nonlinear Ocean Structures,
  WORLD SCIENTIFIC, 1998.
http://dx.doi.org/10.1142/3717 doi:10.1142/3717.

monin2007b
A. S. Monin, A. M. Yaglom, Statistical Fluid Mechanics, Volume II:
  Mechanics of Turbulence, Dover, 2007.

pope1994
S. B. Pope, Lagrangian PDF methods for turbulent flows, Annu. Rev. Fluid
  Mech. 26 (1) (1994) 23–63.
http://dx.doi.org/10.1146/annurev.fl.26.010194.000323
  doi:10.1146/annurev.fl.26.010194.000323.

frisch1995
U. Frisch, Turbulence: The legacy of A. N. Kolmogorov, Cambridge University
  Press, 1995.

venturi2018
D. Venturi, The numerical approximation of nonlinear functionals and functional
  differential equations, Phys. Rep. 732 (2018) 1–102.
http://dx.doi.org/10.1016/j.physrep.2017.12.003
  doi:10.1016/j.physrep.2017.12.003.

cowen2017
L. Cowen, T. Ideker, B. J. Raphael, R. Sharan, Network propagation: A
  universal amplifier of genetic associations, Nat Rev Genet 18 (9) (2017)
  551–562.
http://dx.doi.org/10.1038/nrg.2017.38
  doi:10.1038/nrg.2017.38.

justin2002
J. Zinn-Justin, Quantum Field Theory and Critical Phenomena, 4th Edition,
  Oxford University Press, 2002.
http://dx.doi.org/10.1093/acprof:oso/9780198509233.001.0001
  doi:10.1093/acprof:oso/9780198509233.001.0001.

amit2005
D. J. Amit, V. Martin-Mayor, Field Theory, the Renormalization Group, and
  Critical Phenomena, WORLD SCIENTIFIC, 2005.
http://dx.doi.org/10.1142/5715 doi:10.1142/5715.

risken1989
H. Risken, The Fokker-Planck Equation, 2nd Edition, Springer Berlin
  Heidelberg, 1989.
http://dx.doi.org/10.1007/978-3-642-61544-3
  doi:10.1007/978-3-642-61544-3.

tartakovsky2015
D. M. Tartakovsky, P. A. Gremaud, Method of distributions for uncertainty
  quantification, in: R. Ghanem, D. Higdon, H. Owhadi (Eds.), Handbook of
  Uncertainty Quantification, Springer International Publishing, Switzerland,
  2017, pp. 763–783.
http://dx.doi.org/10.1007/978-3-319-12385-1_27
  doi:10.1007/978-3-319-12385-1_27.

dimarco2014
G. Dimarco, L. Pareschi, Numerical methods for kinetic equations, Acta Numerica
  23 (2014) 369–520.
http://dx.doi.org/10.1017/s0962492914000063
  doi:10.1017/s0962492914000063.

rjasanow2004
S. Rjasanow, W. Wagner, Stochastic numerics for the Boltzmann equation,
  Springer, 2004.

Babovsky1986
H. Babovsky, H. Neunzert, On a simulation scheme for the Boltzmann equation,
  Math. Meth. Appl. Sci. 8 (1) (1986) 223–233.
http://dx.doi.org/10.1002/mma.1670080114
  doi:10.1002/mma.1670080114.

ansumali2011
S. Ansumali, Mean-field model beyond Boltzmann-enskog picture for dense
  gases, Commun. Commut. Phys. 9 (05) (2011) 1106–1116.
http://dx.doi.org/10.4208/cicp.301009.240910s
  doi:10.4208/cicp.301009.240910s.

ramanathan2009
S. Ramanathan, D. L. Koch, An efficient direct simulation Monte Carlo
  method for low mach number noncontinuum gas flows based on the
  bhatnagar–gross–krook model, Phys. Fluids 21 (3) (2009) 033103.
http://dx.doi.org/10.1063/1.3081562
  doi:10.1063/1.3081562.

peraud2014
J.-P. M. Peraud, C. D. Landon, N. G. Hadjiconstantinou, Monte Carlo methods
  for solving the Boltzmann transport equation, Annual Rev Heat Transfer
  17 (N/A) (2014) 205–265.
http://dx.doi.org/10.1615/annualrevheattransfer.2014007381
  doi:10.1615/annualrevheattransfer.2014007381.

rogier1994
F. Rogier, J. Schneider, A direct method for solving the Boltzmann equation,
  Transp. Theory Stat. Phys. 23 (1-3) (1994) 313–338.
http://dx.doi.org/10.1080/00411459408203868
  doi:10.1080/00411459408203868.

cho2016
H. Cho, D. Venturi, G. Karniadakis, Numerical methods for high-dimensional
  probability density function equations, J. Comput. Phys. 305 (2016) 817–837.
http://dx.doi.org/10.1016/j.jcp.2015.10.030
  doi:10.1016/j.jcp.2015.10.030.

zhang2017
Z. Zhang, G. E. Karniadakis, Numerical Methods for Stochastic Partial
  Differential Equations with White Noise, Springer International Publishing,
  2017.
http://dx.doi.org/10.1007/978-3-319-57511-7
  doi:10.1007/978-3-319-57511-7.

weinan2017
W. E, J. Han, A. Jentzen, Deep learning-based numerical methods for
  high-dimensional parabolic partial differential equations and backward
  stochastic differential equations, Commun. Math. Stat. 5 (4) (2017) 349–380.
http://dx.doi.org/10.1007/s40304-017-0117-6
  doi:10.1007/s40304-017-0117-6.

Bachmayr2016
M. Bachmayr, R. Schneider, A. Uschmajew, Tensor networks and hierarchical
  tensors for the solution of high-dimensional partial differential equations,
  Found Comput Math 16 (6) (2016) 1423–1472.
http://dx.doi.org/10.1007/s10208-016-9317-9
  doi:10.1007/s10208-016-9317-9.

hesthaven2007
J. S. Hesthaven, S. Gottlieb, D. Gottlieb, Spectral Methods for Time-Dependent
  Problems, Vol. 21, Cambridge University Press, 2007.
http://dx.doi.org/10.1017/cbo9780511618352
  doi:10.1017/cbo9780511618352.

beylkin2009
G. Beylkin, J. Garcke, M. J. Mohlenkamp, Multivariate regression and machine
  learning with sums of separable functions, SIAM J. Sci. Comput. 31 (3) (2009)
  1840–1857.
http://dx.doi.org/10.1137/070710524
  doi:10.1137/070710524.

acar2011
E. Acar, D. M. Dunlavy, T. G. Kolda, A scalable optimization approach for
  fitting canonical tensor decompositions, J. Chemometrics 25 (2) (2011)
  67–86.
http://dx.doi.org/10.1002/cem.1335
  doi:10.1002/cem.1335.

espig2012
M. Espig, W. Hackbusch, A regularized newton method for the efficient
  approximation of tensors represented in the canonical tensor format, Numer.
  Math. 122 (3) (2012) 489–525.
http://dx.doi.org/10.1007/s00211-012-0465-9
  doi:10.1007/s00211-012-0465-9.

karlsson2016
L. Karlsson, D. Kressner, A. Uschmajew, Parallel algorithms for tensor
  completion in the CP format, Parallel Comput. 57 (2016) 222–234.
http://dx.doi.org/10.1016/j.parco.2015.10.002
  doi:10.1016/j.parco.2015.10.002.

doostan2013
A. Doostan, A. Validi, G. Iaccarino, Non-intrusive low-rank separated
  approximation of high-dimensional stochastic models, Comput. Methods Appl.
  Mech. Eng. 263 (2013) 42–55.
http://dx.doi.org/10.1016/j.cma.2013.04.003
  doi:10.1016/j.cma.2013.04.003.

reynolds2016
M. J. Reynolds, A. Doostan, G. Beylkin, Randomized alternating least squares
  for canonical tensor decompositions: Application to a PDE with random
  data, SIAM J. Sci. Comput. 38 (5) (2016) A2634–A2664.
http://dx.doi.org/10.1137/15m1042802
  doi:10.1137/15m1042802.

battaglino2017
C. Battaglino, G. Ballard, T. G. Kolda, A practical randomized CP tensor
  decomposition (2017).
http://arxiv.org/abs/1701.06600 arXiv:1701.06600.

uschmajew2012
A. Uschmajew, Local convergence of the alternating least squares algorithm for
  canonical tensor approximation, SIAM J. Matrix Anal. & Appl. 33 (2) (2012)
  639–652.
http://dx.doi.org/10.1137/110843587
  doi:10.1137/110843587.

espig2015
M. Espig, W. Hackbusch, A. Khachatryan, On the convergence of alternating least
  squares optimisation in tensor format representations (2015).
http://arxiv.org/abs/1506.00062 arXiv:1506.00062.

bezdek2003
J. C. Bezdek, R. J. Hathaway, Convergence of alternating optimization, Neural
  Parallel Sci. Comput. 11 (2003) 351–368.

rohwedder2013
T. Rohwedder, A. Uschmajew, On local convergence of alternating schemes for
  optimization of convex problems in the tensor train format, SIAM J. Numer.
  Anal. 51 (2) (2013) 1134–1162.
http://dx.doi.org/10.1137/110857520
  doi:10.1137/110857520.

ortega1970
J. M. Ortega, W. C. Rheinboldt, Iterative Solution of Nonlinear Equations in
  Several Variables, Society for Industrial and Applied Mathematics, 2000.
http://dx.doi.org/10.1137/1.9780898719468
  doi:10.1137/1.9780898719468.

kaya2018
O. Kaya, B. Uçar, Parallel Candecomp/Parafac decomposition of sparse tensors
  using dimension trees, SIAM J. Sci. Comput. 40 (1) (2018) C99–C130.
http://dx.doi.org/10.1137/16m1102744
  doi:10.1137/16m1102744.

Hackbusch2009
W. Hackbusch, S. Kühn, A new scheme for the tensor representation, J Fourier
  Anal Appl 15 (5) (2009) 706–722.
http://dx.doi.org/10.1007/s00041-009-9094-9
  doi:10.1007/s00041-009-9094-9.

hackbusch2012
W. Hackbusch, Tensor Spaces and Numerical Tensor Calculus, Springer Berlin
  Heidelberg, 2012.
http://dx.doi.org/10.1007/978-3-642-28027-6
  doi:10.1007/978-3-642-28027-6.

kolda2009
T. G. Kolda, B. W. Bader, Tensor decompositions and applications, SIAM Rev.
  51 (3) (2009) 455–500.
http://dx.doi.org/10.1137/07070111x
  doi:10.1137/07070111x.

grasedyck2010
L. Grasedyck, Hierarchical singular value decomposition of tensors, SIAM J.
  Matrix Anal. & Appl. 31 (4) (2010) 2029–2054.
http://dx.doi.org/10.1137/090764189
  doi:10.1137/090764189.

lathauwer2000
L. De Lathauwer, B. De Moor, J. Vandewalle, A multilinear singular value
  decomposition, SIAM J. Matrix Anal. & Appl. 21 (4) (2000) 1253–1278.
http://dx.doi.org/10.1137/s0895479896305696
  doi:10.1137/s0895479896305696.

venturi2011
D. Venturi, A fully symmetric nonlinear biorthogonal decomposition theory for
  random fields, Physica D 240 (4-5) (2011) 415–425.
http://dx.doi.org/10.1016/j.physd.2010.10.005
  doi:10.1016/j.physd.2010.10.005.

peres1995
A. Peres, Higher order Schmidt decompositions, Phys. Lett. A 202 (1) (1995)
  16–17.
http://dx.doi.org/10.1016/0375-9601(95)00315-t
  doi:10.1016/0375-9601(95)00315-t.

silva2008
V. de Silva, L.-H. Lim, Tensor rank and the ill-posedness of the best low-rank
  approximation problem, SIAM J. Matrix Anal. & Appl. 30 (3) (2008)
  1084–1127.
http://dx.doi.org/10.1137/06066518x
  doi:10.1137/06066518x.

hillar2013
C. J. Hillar, L.-H. Lim, Most tensor problems are NP-hard, JACM 60 (6) (2013)
  1–39.
http://dx.doi.org/10.1145/2512329
  doi:10.1145/2512329.

vannieuwenhoven2014
N. Vannieuwenhoven, J. Nicaise, R. Vandebril, K. Meerbergen, On generic
  nonexistence of the Schmidt–eckart–young decomposition for complex
  tensors, SIAM J. Matrix Anal. & Appl. 35 (3) (2014) 886–903.
http://dx.doi.org/10.1137/130926171
  doi:10.1137/130926171.

Kressner2014
D. Kressner, C. Tobler, Algorithm 941, ACM Trans. Math. Softw. 40 (3) (2014)
  1–22.
http://dx.doi.org/10.1145/2538688
  doi:10.1145/2538688.

DaSilva2015
C. Da Silva, F. J. Herrmann, Optimization on the hierarchical tucker manifold
  – applications to tensor completion, Linear Algebra and its Applications 481
  (2015) 131–173.
http://dx.doi.org/10.1016/j.laa.2015.04.015
  doi:10.1016/j.laa.2015.04.015.

oseledets2011
I. V. Oseledets, Tensor-train decomposition, SIAM J. Sci. Comput. 33 (5) (2011)
  2295–2317.
http://dx.doi.org/10.1137/090752286
  doi:10.1137/090752286.

grasedyck2015
L. Grasedyck, R. Kriemann, C. Löbbert, A. Nägel, G. Wittum, K. Xylouris,
  Parallel tensor sampling in the hierarchical tucker format, Comput. Visual
  Sci. 17 (2) (2015) 67–78.
http://dx.doi.org/10.1007/s00791-015-0247-x
  doi:10.1007/s00791-015-0247-x.

Nouy2017
A. Nouy, Higher-order principal component analysis for the approximation of
  tensors in tree-based low-rank formats (2017).
http://arxiv.org/abs/1705.00880 arXiv:1705.00880.

etter2016
S. Etter, Parallel ALS algorithm for solving linear systems in the
  hierarchical tucker representation, SIAM J. Sci. Comput. 38 (4) (2016)
  A2585–A2609.
http://dx.doi.org/10.1137/15m1038852
  doi:10.1137/15m1038852.

grasedyck2017
L. Grasedyck, C. Löbbert, Distributed Hierarchical SVD in the
  Hierarchical Tucker Format (2017).
http://arxiv.org/abs/1708.03340 arXiv:1708.03340.

Phan2013b
A.-H. Phan, P. Tichavsky, A. Cichocki, CANDECOMP/PARAFAC decomposition of
  high-order tensors through tensor reshaping, IEEE Trans. Signal Process.
  61 (19) (2013) 4847–4860.
http://dx.doi.org/10.1109/tsp.2013.2269046
  doi:10.1109/tsp.2013.2269046.

Zhou2013
G. Zhou, A. Cichocki, S. Xie, Accelerated canonical polyadic decomposition
  using mode reduction, IEEE Trans. Neural Netw. Learning Syst. 24 (12) (2013)
  2051–2062.
http://dx.doi.org/10.1109/tnnls.2013.2271507
  doi:10.1109/tnnls.2013.2271507.

Hatch2012
D. Hatch, D. del Castillo-Negrete, P. Terry, Analysis and compression of
  six-dimensional gyrokinetic datasets using higher order singular value
  decomposition, J. Comput. Phys. 231 (11) (2012) 4234–4256.
http://dx.doi.org/10.1016/j.jcp.2012.02.007
  doi:10.1016/j.jcp.2012.02.007.

Kormann2015
K. Kormann, A semi-Lagrangian vlasov solver in tensor train format, SIAM J.
  Sci. Comput. 37 (4) (2015) B613–B632.
http://dx.doi.org/10.1137/140971270
  doi:10.1137/140971270.

Dolgov2014
S. Dolgov, A. Smirnov, E. Tyrtyshnikov, Low-rank approximation in the numerical
  modeling of the Farley–buneman instability in ionospheric plasma, J.
  Comput. Phys. 263 (2014) 268–282.
http://dx.doi.org/10.1016/j.jcp.2014.01.029
  doi:10.1016/j.jcp.2014.01.029.

quarteroni2007
A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics, 2nd Edition,
  Springer New York, 2007.
http://dx.doi.org/10.1007/978-0-387-22750-4
  doi:10.1007/978-0-387-22750-4.

rhee2001
H.-K. Rhee, R. Aris, N. R. Amundson, First-order partial differential
  equations, volume 1: Theory and applications of single equations, Dover,
  2001.

dimarco2018
G. Dimarco, R. Loubère, J. Narski, T. Rey, An efficient numerical method for
  solving the Boltzmann equation in multidimensions, J. Comput. Phys. 353
  (2018) 46–81.
http://dx.doi.org/10.1016/j.jcp.2017.10.010
  doi:10.1016/j.jcp.2017.10.010.

desvillettes1996
L. Desvillettes, S. Mischler, About the splitting algorithm for Boltzmann and
  B.G.K. equations, Math. Models Methods Appl. Sci. 06 (08) (1996)
  1079–1101.
http://dx.doi.org/10.1142/s0218202596000444
  doi:10.1142/s0218202596000444.

shuleshko1959
P. Shuleshko, A new method of solving boundary-value problems of mathematical
  physics, Aust. J. Appl. Sci 10 (1959) 1–7.

snyder1964
L. J. Snyder, T. W. Spriggs, W. E. Stewart, Solution of the equations of change
  by galerkin's method, AIChE J. 10 (4) (1964) 535–540.
http://dx.doi.org/10.1002/aic.690100423
  doi:10.1002/aic.690100423.

zinn1968
B. Zinn, E. Powell, Application of the Galerkin method in the solution of
  combustion-instability problems, in: XlXth International Astronautical
  Congress, Vol. 3, 1968, pp. 59–73.

comon2009
P. Comon, X. Luciani, A. L. F. de Almeida, Tensor decompositions, alternating
  least squares and other tales, J. Chemometrics 23 (7-8) (2009) 393–405.
http://dx.doi.org/10.1002/cem.1236
  doi:10.1002/cem.1236.

cercignani1997
C. Cercignani, V. I. Gerasimenko, D. Y. Petrina, Many-Particle Dynamics and
  Kinetic Equations, 1st Edition, Springer Netherlands, 1997.
http://dx.doi.org/10.1007/978-94-011-5558-8
  doi:10.1007/978-94-011-5558-8.

villani2002
C. Villani, A review of mathematical topics in collisional kinetic theory,
  in: S. Friedlander, D. Serre (Eds.), Handbook of mathematical fluid
  mechanics, Vol. 1, North-Holland, 2002, pp. 71–305.

cercignani1994
C. Cercignani, R. Illner, M. Pulvirenti, The Mathematical Theory of Dilute
  Gases, Springer New York, 1994.
http://dx.doi.org/10.1007/978-1-4419-8524-8
  doi:10.1007/978-1-4419-8524-8.

Nishida1978
T. Nishida, Fluid dynamical limit of the nonlinear Boltzmann equation to the
  level of the compressible Euler equation, Commun.Math. Phys. 61 (2) (1978)
  119–148.
http://dx.doi.org/10.1007/bf01609490
  doi:10.1007/bf01609490.

caflisch1980
R. E. Caflisch, The fluid dynamic limit of the nonlinear boltzmann equation,
  Comm. Pure Appl. Math. 33 (5) (1980) 651–666.
http://dx.doi.org/10.1002/cpa.3160330506
  doi:10.1002/cpa.3160330506.

Struchtrup2005
H. Struchtrup, Macroscopic Transport Equations for Rarefied Gas Flows, Springer
  Berlin Heidelberg, 2005.
http://dx.doi.org/10.1007/3-540-32386-4
  doi:10.1007/3-540-32386-4.

Levermore1996
C. D. Levermore, Moment closure hierarchies for kinetic theories, J Stat Phys
  83 (5-6) (1996) 1021–1065.
http://dx.doi.org/10.1007/bf02179552
  doi:10.1007/bf02179552.

mueller1998
I. Müller, T. Ruggeri, Rational Extended Thermodynamics, Springer New York,
  1998.
http://dx.doi.org/10.1007/978-1-4612-2210-1
  doi:10.1007/978-1-4612-2210-1.

bhatnagar1954
P. L. Bhatnagar, E. P. Gross, M. Krook, A model for collision processes in
  gases. i. small amplitude processes in charged and neutral one-component
  systems, Phys. Rev. 94 (3) (1954) 511–525.
http://dx.doi.org/10.1103/physrev.94.511
  doi:10.1103/physrev.94.511.

Mieussens2000
L. Mieussens, Discrete-velocity models and numerical schemes for the
  Boltzmann-BGK equation in plane and axisymmetric geometries, J. Comput.
  Phys. 162 (2) (2000) 429–466.
http://dx.doi.org/10.1006/jcph.2000.6548
  doi:10.1006/jcph.2000.6548.

nassios2013
J. Nassios, J. E. Sader, High frequency oscillatory flows in a
  slightly~rarefied gas according to the
  Boltzmann–BGK~equation, J. Fluid Mech. 729 (2013) 1–46.
http://dx.doi.org/10.1017/jfm.2013.281
  doi:10.1017/jfm.2013.281.

Andries2000
P. Andries, P. Le Tallec, J.-P. Perlat, B. Perthame, The Gaussian-BGK model
  of Boltzmann equation with small prandtl number, Eur. J. Mech. B. Fluids
  19 (6) (2000) 813–830.
http://dx.doi.org/10.1016/s0997-7546(00)01103-1
  doi:10.1016/s0997-7546(00)01103-1.

borghini2016
N. Borghini, Topics in Nonequilibrium Physics (Sep. 2016).

barichello2003
L. Barichello, C. Siewert, Some comments on modeling the linearized Boltzmann
  equation, J. Quant. Spectrosc. Radiat. Transfer 77 (1) (2003) 43–59.
http://dx.doi.org/10.1016/s0022-4073(02)00074-2
  doi:10.1016/s0022-4073(02)00074-2.

pekeris1957
C. L. Pekeris, Z. Alterman, Solution of the Boltzmann-Hilbert integral
  equation ii. the coefficients of viscosity and heat conduction, Proceedings
  of the National Academy of Sciences 43 (11) (1957) 998–1007.
http://dx.doi.org/10.1073/pnas.43.11.998
  doi:10.1073/pnas.43.11.998.

loyalka1967
S. K. Loyalka, Model dependence of the slip coefficient, Phys. Fluids 10 (8)
  (1967) 1833.
http://dx.doi.org/10.1063/1.1762366
  doi:10.1063/1.1762366.

huang2013
H. Huang, J. M. Dennis, L. Wang, P. Chen, A scalable parallel LSQR algorithm
  for solving large-scale linear system for tomographic problems: A case
  study in seismic tomography, Procedia Comput. Sci. 18 (2013) 581–590.
http://dx.doi.org/10.1016/j.procs.2013.05.222
  doi:10.1016/j.procs.2013.05.222.

paige1982
C. C. Paige, M. A. Saunders, LSQR: An algorithm for sparse linear equations
  and sparse least squares, ACM Trans. Math. Softw. 8 (1) (1982) 43–71.
http://dx.doi.org/10.1145/355984.355989
  doi:10.1145/355984.355989.

barrett1994
R. Barrett, M. Berry, T. F. Chan, J. Demmel, J. Donato, J. Dongarra,
  V. Eijkhout, R. Pozo, C. Romine, H. van der Vorst, Templates for the Solution
  of Linear Systems: Building Blocks for Iterative Methods, Society for
  Industrial and Applied Mathematics, 1994.
http://dx.doi.org/10.1137/1.9781611971538
  doi:10.1137/1.9781611971538.






