In [59]:
from IPython.html.services.config import ConfigManager
from IPython.utils.path import locate_profile
cm = ConfigManager(profile_dir=locate_profile(get_ipython().profile))
cm.update('livereveal', {
'theme': 'sky',
'transition': 'zoom',
'start_slideshow_at': 'selected',
})
%load_ext tikzmagic
The numerical methods for PDEs are well-established:
Typically, the "order" of the numerical scheme is considered, i.e.
$$\Vert u - u_h \Vert = \mathcal{O}(h^p),$$where $h$ is the grid size
We write down the weak formulation for the equation,
$$-\int \nabla u \nabla v dx = \int f v dx, \quad v \in W^{2}_1,$$and then use the Galerkin method, replacing $u$ by $\sum_i u_i \varphi_i$ and testing the weak equation with $v = \phi_j$.
You can do it automatically using the FeNics or Firedrake software.
In [76]:
from firedrake import *
n = 128
mesh = UnitSquareMesh(n, n)
V = FunctionSpace(mesh, "CG", 2)
bc0 = DirichletBC(V, Expression(("0")), 1)
bc1 = DirichletBC(V, Expression(("0")), 2)
bc2 = DirichletBC(V, Expression(("0")), 3)
#bc3 = DirichletBC(V, Expression(("0")), 4)
bcs = [bc0, bc1, bc2]
u = TrialFunction(V)
v = TestFunction(V)
f = Function(V)
f.interpolate(Expression("1"))
#a = ( dot(grad(v), grad(u))) * dx
#frm = assemble(a, bcs=bcs)
L = f * v * dx
L1 = assemble(L)
a = (dot(grad(v), grad(u))) * dx
frm = assemble(a, bcs=bcs)
sol = Function(V)
solve(frm, sol, L1, solver_parameters={'ksp_type': 'cg'})
sol.vector().array()
#sol = Function(V)
#solve(a, sol, L1, solver_parameters={'ksp_type': 'cg'})
Out[76]:
In [67]:
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns
u0 = sol.vector().array()
p1 = mesh.coordinates.vector().array().reshape((-1, 2))
plt.tricontourf(p1[:, 0], p1[:, 1], u0)
What are the challenges, if FEM (and also other techniques, like Finite Volume, Finite Difference methods) are so well-developed?
There are several:
There are different ways to adapt the solution to get better accuracy for the same number of degrees of freedom.
Hp-refinement is the most efficient method, but it is technically very non-trivial to implement,
and requires the knowledge of where the singularity is, what is the order of the singularity of the solution,
where $k = k(\frac{x}{\delta})$ -- periodic case.
Is to expand the solution in "series" of small parameter $\delta$.
The "zero"-order approximation is to replace $k$ by the averaged coefficient over the unit cell:
$$ \nabla \cdot (k_0 \nabla u_0) = f,$$and $k_0$ is the effective coefficient.
Next terms include expansion of the solution as
$$ u(x, \frac{x}{\delta}) = u_0 + \delta \sigma(\delta^{-1} \rho(x)) u_1\left(x, \frac{x}{\delta}\right), $$where $\sigma(0) = 0, 0 \leq \sigma(t) \leq 1$ for $0 \leq t \leq 1$ and $\sigma(t) = 1$ for all $t \geq 1$.
$\rho$ is the distance to the boundary.
For many problems, the asymptotic expansions are written. They include boundary layers, quasi-periodic coefficients, vector problems.
Computational complexity is solving several unit cells problems (with different boundary conditions).
Main problem: a lot of analytic work, not 100% flexible (i.e., you want to make a defect in one point)
In [70]:
#Generate sparse tridiag matrix for div K grad
import scipy.sparse
import numpy as np
def Kfun(x, eps):
return (1 + x) * 2/3 * (1 + np.cos(2 * np.pi * (x/eps)) ** 2)
def gen_full_matrix(d, eps, Kfun, a=0.0, b=1.0):
N = 2 ** d
t = np.arange(N)
e = np.ones(N)
h = (b - a)/(N + 1)
cf = (t + 0.5 * e) * h
am = Kfun(cf, eps)
cf_plus = (t + 1.5 * e) * h
aa = (np.arange(N+1) + 0.5 * np.ones(N+1)) * h
aa = Kfun(aa, eps)
ap = Kfun(cf_plus, eps)
dg = am + ap
mat = scipy.sparse.spdiags([-am, -ap, dg], [1, -1, 0], N, N)
mat = mat/h**2
return aa, am, ap, mat
d = 15
eps = 1e-3
mat = gen_full_matrix(d, eps, Kfun)[-1]
rhs = np.ones(mat.shape[0])
sol = scipy.sparse.linalg.spsolve(mat, rhs)
fig, ax = plt.subplots(1, 2, figsize=(10, 4))
ax[0].plot(sol)
ax[1].plot(np.diff(sol))
ax[0].set_title('Solution plot')
ax[1].set_title('Derivative of the solution')
fig.tight_layout()
The idea of QTT is very simple (and has been around for many many years).
Suppose we have a one-dimensional function, that is sampled on a uniform grid with $2^d$ points.
$$v_k = f(x_k), \quad k = 1, \ldots, 2^d$$Then we reshape this vector into $2 \times 2 \times \ldots \times 2$ $d$-dimensional tensor:
$$V(i_1, \ldots, i_d) = v(i), \quad i = i_1+ 2 i_1 + \ldots 2^{d-1} i_d.$$Given a tensor $V(i_1, \ldots, i_d)$ we say that it is in the canonical format, if it is the sum of rank-1 tensors:
$$V(i_1, \ldots, i_d) = \sum_{\alpha=1}^r V_1(i_1, \alpha) \ldots V_d(i_d, \alpha).$$For $d=2$ it reduces to the low-rank approximation of a matrix, thus for many years it was considered as a multidimensional generalization of the singular value decomposition.
Canonical decomposition has problems:
However, there exists another multidimensional generalization of the SVD, which may have larger number of parameters, but it is much easier to compute!
Tensor-train decomposition has the form
$$A(i_1, \ldots, i_d) = G_1(i_1) \ldots G_d(i_d),$$where $G_k(i_k)$ is $r_{k-1} \times r_k$ for any fixed $i_k$, and $r_0 = r_d = 1$.
It is also known for many years as matrix product states in solid state physics (and used to approximate wavefunctions).
QTT decomposition = Quantization + TT
I.e. we take a function, sample it on a uniform grid, replace it by a $d$-dimensional tensor, approximate it by the TT-format.
We get $\mathcal{O}(d r^2)$ instead of $2^d$ parameters, but this is a non-linear representation.
Main question: what kind of functions can we represent in such way?
In two dimensions, if the domain is a square, the situation is the same:
Discretize on a tensor-product mesh, quantize each one-dimensional index and go.
But what to do for non-square domains?
One idea is to use characteristic functions, but it does not work.
The solution, proposed by Vladimir Kazeev and Christoph Schwab is astonishingly simple and efficient:
Introduce logically rectangular grids!
Kazeev and Schwab in the paper Quantized tensor-structured finite elements for second-order elliptic PDEs in two dimensions
have proven the following theoretical result, which I very loosely formulate as follows.
For a second-order PDE with analytic forcing, and for a given number of parameters $N$ in the QTT representation, there exists a QTT-FEM approximant such that
$$\Vert u - u_{QTT, h} \Vert_{\mathbb{H}_1} = \mathcal{O}(e^{-bN^{\alpha}}),$$where $\alpha=\frac{1}{5}$ in their theoretical analysis (in experiments it is much better, $\alpha \approx \frac{1}{2}$).
Idea of the proof is to use known results on the hp-method and on the regularity of the solutions in polygonal domains.
Singularities have very specific form (point singularities), which can be approximated in the QTT-format.
The remainder is the analytic function, which can be approximated by the Fourier series (and they have small QTT-rank).
We can take a low-order FEM on a very fine virtual mesh.
Formulate a huge linear system
$$Ax = f,$$and look for the solution in the TT-format.
The general concept is very simple: we reformulate the solution as a minimization of energy:
$$(Ax, x) - 2(f, x) \rightarrow \min,$$and minimize over the set of tensors with bounded ranks.
This is non-convex minimization procedure with low-rank constraints, but with much fewer number of parameters, and we have developed a black-box software for solving such problems.
TT-Toolbox is our main software packages, available:
In [71]:
import tt
from tt.amen import amen_solve
d = 8
a = tt.qlaplace_dd([d, d, d])
rhs = tt.ones(2, 3*d)
sol = amen_solve(a, rhs, tt.rand(rhs.n), 1e-6) #Black-box solver
mem = int((sol.erank ** 2) * sum(sol.n))
print('Full memory: {0: d}, QTT-memory: {1: d}'.format(np.prod(sol.n), mem))
Now let us get back to the multiscale problems.
There are two important examples of multiscale problems:
The concept is the same: introduce a virtual mesh that resolves the finest scale, approximate it.
We have recently proved (Kazeev, Schwab, Oseledets, Rakhuba) that the same result hold for 1D multiscale problems!
Given the equation
$$(k u')' = f, $$where $K = K\left(x, \frac{x}{\delta}\right)$
the solution can be approximated by a QTT-interpolant with exponential convergence with respect to the number of parameters in the QTT-format. The number of parameters required to reach the accuracy is also logarithmic in the multiscale parameter $\delta$, i.e.
$$N_l \sim \log^{\alpha} \epsilon \log^{\beta} \delta.$$
In [1]:
import numpy as np
import tt
import tt.amen
from tt.amen import amen_solve
import time
from tt.multifuncrs2 import multifuncrs2
d = 50
N = 2 ** d
x = tt.ones(2, d)
h = 1.0/N
B = (tt.Toeplitz(x, kind='U') + tt.eye(2, d))#.full()
B = B * h
alpha = -0.5
st = time.time()
tau = 1e-8
t = tt.xfun(2, d)
x = t * h
eps = 1e-4
Kf = lambda x: Kfun(x, eps)
Kf = multifuncrs2([x], Kf, tau, verb=0)
Kf = Kf.round(tau)
Kmat = tt.diag(Kf)
ed = tt.eye(2, d)
#idd = tt.matrix(np.array([[0, -1.0], [0, 0]]))
e00 = tt.matrix(np.array([[1.0, 0, 0], [0, 0, 0], [0, 0, 0]]))
e01 = tt.matrix(np.array([[0.0, -1.0, 0], [0, 0, 0], [0, 0, 0]]))
e11 = tt.matrix(np.array([[0, 0, 0], [0, 1.0, 0], [0, 0, 0]]))
e12 = tt.matrix(np.array([[0, 0, 0], [0, 0.0, -1.0], [0, 0, 0]]))
e22 = tt.matrix(np.array([[0, 0, 0], [0, 0.0, 0.0], [0, 0, 1.0]]))
big_mat = tt.kron(ed, e00) + tt.kron(B, e01) + tt.kron(Kmat, e11) + tt.kron(ed, e12) + tt.kron(ed, e22)
f = tt.ones(2, d)
alpha = -0.5
rhs1 = tt.matvec(B, f) + alpha * tt.ones(2, d)
ev = tt.tensor(np.array([0.0, 0.0, 1.0]))
rhs = tt.kron(rhs1, ev)
sol = amen_solve(big_mat, rhs, rhs, 1e-6)
st = time.time() - st
print 'Total computing time:', st
#sol_full = sol.full().flatten("F")
#u1 = sol_full[:N]
#w1 = sol_full[N:2*N]
#fig, ax = plt.subplots(1, 3, figsize=(12, 4))
#x0 = x.full().flatten("F")
#ax[0].plot(x0, u1)
#ax[0].set_title('Solution')
#ax[1].plot(x0, w1)
#ax[1].set_title('Derivative')
#ax[2].plot(w1[:60])
#ax[2].set_title('Zooming in the derivative')
#fig.text(0.5, -0.1, 'Eps = {0}, d = {1}'.format(eps, d), size=16)
#fig.tight_layout()
#ax[1]
Suppose we have a 1D Laplacian
$$\frac{u_{i-1} - 2 u_i + u_{i+1}}{h^2},$$and the error in this is $\frac{\epsilon}{h^2} + h^2$, i.e. we can not go beyond $h \sim \epsilon^{1/4}.$
For a non-QTT-world, there is no problem: $\epsilon = 10^{-16}, \quad h \sim 10^{-4}$.
In our case, we need grids $2^{60}-2^{80}$, so this is a problem.
We need new discretization schemes (and I already used one in the example above).
Our main goal is to establish the technology at all levels and show its efficiency in real applications. We now have the theoretical ground.
In [58]:
from IPython.core.display import HTML
def css_styling():
styles = open("custom.css", "r").read()
return HTML(styles)
css_styling()
Out[58]: