In [1]:
from traitlets.config.manager import BaseJSONConfigManager
#path = "/home/damian/miniconda3/envs/rise_latest/etc/jupyter/nbconfig"
cm = BaseJSONConfigManager()
cm.update('livereveal', {
'theme': 'sky',
'transition': 'zoom',
'start_slideshow_at': 'selected',
'scroll': True
})
Out[1]:
In the forward propagation of uncertainty, we have a known model F for a system of interest.
We model its inputs $X$ as a random variable and wish to understand the output random variable
$$Y = F(X)$$(also denoted $Y \vert X$) and reads $Y$ given $X$.
We want to infer statistical properties of $Y$ given statistical properties of $X$.
In the simplest case, this correponds to the computation of high-dimensional integrals.
$$E Y = \int F(X) \rho(X) dX.$$There are not too many efficient methods for working with multivariate functions:
Many successful methods are based on the idea of separation of variables
Function is called separable, if $$ f(x_1, \ldots, x_d) = f_1(x_1) f_2(x_2) \ldots f_d(x_d).$$
Not many functions can be represented in this form
(or even approximated, it means independence of variables in some sense.)
A more general and efficient class of functions is the sum of separable functions (known also proper generalized decomposition).
$$f(x_1,\ldots, x_d) \approx \sum_{\alpha=1}^r f^{(\alpha)}_1(x_1) \ldots f^{(\alpha)}_d(x_d)$$In the discrete case, representation $$A(i_1, \ldots, i_d) \approx \sum_{\alpha=1}^r U_1(i_1, \alpha) \ldots U_d(i_d, \alpha)$$ is called canonical format, or canonical-polyadic (CP-format).
The number of parameters is $dnr$, so if $r$ is small, a good idea is that it is possible to fit those parameters.
Cross approximation is a heuristic sampling technique for the selection of "good" rows and columns in a low-rank matrix.
The goal is to find rows and column that maximize the volume of the submatrix.
The maximal volume principle (Goreinov, Tyrtyshnikov) states that
$$ \Vert A - A_{skel} \Vert_C \leq (r + 1) \sigma_{r+1},$$and $\Vert \cdot \Vert$ is the maximal in modulus element of a matrix.
You can compute Tucker by means of the SVD:
The simplest construction separates one index at a time and gives the tensor-train (or matrix-product state) decomposition
$$A(i_1, \ldots, i_d) \approx G_1(i_1) \ldots G_d(i_d), $$where $G_k(i_k)$ is $r_{k-1} \times r_k$ matrix for any fixed $r_k$.
Similar to Hidden Markov models
Theorem: there exists a TT-decomposition with TT-ranks
$$r_k = \mathrm{rank} A_k$$
Rounding procedure solves the following problem: given a TT-representation
$$A(i_1, \ldots, i_d) = G_1(i_1) G_2(i_2) \ldots G_d(i_d)$$find another one with smaller memory that approximates it with threshold $\epsilon$.
It can be done in $\mathcal{O}(dnr^3)$ operations using the linear structure of the format:
We can do orthogonalization of the cores!
We have two basic packages:
Many basic and advanced algorithms are implemented, and can be used as building blocks.
In [15]:
import numpy as np
import tt
from tt.cross.rectcross import cross
d = 20
n = 10
x0 = tt.rand(d, n, 3)
h = 1e-2
def myfun1(x):
return np.sum(x, axis=1)
#return np.sqrt((((h * x) ** 2).sum(axis=1)))
x1 = cross(myfun1, x0)
x1 = x1.round(1e-8)
print(x1)
Sometime the function is given only implicitly. How to construct an approximation?
General strategy: formulate the problem as a minimization problem
$$ F(X) \rightarrow \min, \quad X \in \mathcal{M}. $$
In [18]:
import numpy as np
import tt
from tt.amen import amen_solve
d = 8
f = 4
mat = tt.qlaplace_dd([d] * f)
rhs = tt.ones(2, d * f)
sol = amen_solve(mat, rhs, rhs, 1e-6)
print('Error:', (tt.matvec(mat, sol) - rhs).norm()/rhs.norm())
print('Efficient rank:', sol.erank)