概述
在做数值实验时很可能要求实现符合某种协方差的随机域。选取这样的随机向量的一种方式就是Karhunen-Loève expansion。本文主要参考论文介绍基本的概念和简要介绍数值实现的方法。为以后使用做准备。
文章目录
- 1. 随机域
- 2. KL expansion
- 截断KL expansion
- 3. 解KL expansion的数值方法
- Degenerate kernel methods
- Nyström methods
- Projection methods
- (1) Collocation methods
- (2) Galerki methods
1. 随机域
Ω ⊂ R d Omegasubset mathbb{R}^{d} Ω⊂Rd,是一个连续的区域。连续随机域(continuous random field) H ( x , θ ) H(mathbf{x},theta) H(x,θ) 是定义在 Ω Omega Ω上的随机函数, θ ∈ Θ thetainTheta θ∈Θ, ( Θ , F , P ) (Theta,F,P) (Θ,F,P)是一个完备的随机空间。
如果每个 x ∈ Ω mathbf{x}inOmega x∈Ω对应的都是随机变量,那么这个随机域就被成为单变量的(univariate) 或者 实值的(real-valued)。如果每个 x mathbf{x} x对应的是随机向量,那么这个随机域就是多变量的(multivariate)。
d = 1 d=1 d=1,随机域称为一维的(one-dimensional), d > 1 d>1 d>1,随机域成为多维的(multidimentional)。
如果 ∀ ( x 1 , … , x n ) ∈ Ω forall (mathbf{x}_1, dots,mathbf{x}_n)inOmega ∀(x1,…,xn)∈Ω, ∀ n forall n ∀n,有 ( H ( x 1 , θ ) , … , H ( x n , θ ) ) (H(mathbf{x}_1,theta),dots,H(mathbf{x}_n,theta)) (H(x1,θ),…,H(xn,θ))是联合高斯分布的,则称随机域是高斯的(Gaussian)。
均值函数: μ : Ω → R mu:Omegato mathbb{R} μ:Ω→R。自相关函数: C o v : Ω × Ω → R mathsf{Cov}:OmegatimesOmegato mathbb{R} Cov:Ω×Ω→R。标准差函数: σ : Ω → R sigma:Omegato mathbb{R} σ:Ω→R。自相关系数函数: ρ : Ω × Ω → [ − 1 , 1 ] rho:OmegatimesOmegato [-1,1] ρ:Ω×Ω→[−1,1]。
C o v ( x , x ′ ) = σ ( x ) σ ( x ′ ) ρ ( x , x ′ ) mathsf{Cov}(mathbf{x},mathbf{x}')=sigma(mathbf{x})sigma(mathbf{x}')rho(mathbf{x},mathbf{x}') Cov(x,x′)=σ(x)σ(x′)ρ(x,x′).
2. KL expansion
Karhunen-Loève expansion是一种以序列形式表达随机域的方式。其建立在自相关函数的谱分解上。 H ( x , θ ) = μ ( x ) + ∑ i = 1 ∞ λ i φ i ( x ) ξ i ( θ ) , H(mathbf{x},theta)=mu(mathbf{x})+sum_{i=1}^{infty}sqrt{lambda_i}varphi_i(mathbf{x})xi_i(theta), H(x,θ)=μ(x)+i=1∑∞λiφi(x)ξi(θ),其中 ∫ Ω C o v ( x , x ′ ) φ i ( x ′ ) d x ′ = λ i φ i ( x ) . (1) int_{Omega}mathsf{Cov}(mathbf{x},mathbf{x}')varphi_{i}(mathbf{x}')mathsf{d}mathbf{x}'=lambda_ivarphi_i(mathbf{x}).tag{1} ∫ΩCov(x,x′)φi(x′)dx′=λiφi(x).(1) φ i varphi_i φi都正规化了,而 ξ i ( θ ) : Θ → R xi_i(theta):Thetato mathbb{R} ξi(θ):Θ→R是标准不相关随机变量。
有效的自相关函数是有界(bounded),对称(symmetric),半正定的(semi-definite)。 λ i ≥ 0 lambda_ige0 λi≥0, C o v ( x , x ′ ) = ∑ i = 1 ∞ λ i φ i ( x ) φ i ( x ′ ) mathsf{Cov}(mathbf{x},mathbf{x}')=sum_{i=1}^{infty}lambda_ivarphi_i(mathbf{x})varphi_i(mathbf{x}') Cov(x,x′)=∑i=1∞λiφi(x)φi(x′), ∫ Ω φ i ( x ) φ ( x ) d x = δ i j int_{Omega}varphi_i(mathbf{x})varphi(mathbf{x})mathsf{d}mathbf{x}=delta_{ij} ∫Ωφi(x)φ(x)dx=δij。 { φ i } {varphi_i} {φi}是 L 2 ( Ω ) mathsf{L}^2(Omega) L2(Ω)的一个完备正交基。
如果 H ( x , θ ) H(mathbf{x},theta) H(x,θ)是高斯的,则 ξ i ( θ ) xi_{i}(theta) ξi(θ)是独立标准高斯随机变量。
截断KL expansion
H
~
(
x
,
θ
)
=
μ
(
x
)
+
∑
i
=
1
M
λ
i
φ
i
(
x
)
ξ
i
(
θ
)
.
tilde{H}(mathbf{x},theta)=mu(mathbf{x})+sum_{i=1}^{M}sqrt{lambda_i}varphi_i(mathbf{x})xi_i(theta).
H~(x,θ)=μ(x)+i=1∑Mλiφi(x)ξi(θ).
截断KL expansion是对原随机域的一种近似。
3. 解KL expansion的数值方法
如果是数值的方式解(1)式,得到的是真实解的近似值 λ ^ i hat{lambda}_i λ^i, φ ^ i hat{varphi}_i φ^i。那么对随机域的近似就成了 H ^ ( x , θ ) = μ ( x ) + ∑ i = 1 M λ ^ i φ ^ i ( x ) ξ ^ i ( θ ) . hat{H}(mathbf{x},theta)=mu(mathbf{x})+sum_{i=1}^{M}sqrt{hat{lambda}_i}hat{varphi}_i(mathbf{x})hat{xi}_{i}(theta). H^(x,θ)=μ(x)+i=1∑Mλ^iφ^i(x)ξ^i(θ).
而 φ ^ i hat{varphi}_i φ^i又是用 { h j : Ω → R } {h_j:Omegato mathbb{R}} {hj:Ω→R}得到的。 φ i ≈ φ ^ i = ∑ j = 1 N d j i h j ( x ) , varphi_iapprox hat{varphi}_i=sum_{j=1}^{N}d_j^{i}h_j(mathbf{x}), φi≈φ^i=j=1∑Ndjihj(x),其中 d i j d_i^{j} dij需要被确定。
解的方法分为以下几类
- Degenerate kernel methods;
- Nyström methods;(这是个德国名字啊,看的这篇论文也是德国人写的,不知道有没有什么关系)
- Projection methods:
(1) Collocation methods;
(2) Galerkin methods.
Degenerate kernel methods
C o v ( x , x ′ ) ≈ K ^ ( x , x ′ ) = ∑ j = 1 N ∑ n = 1 N k j n α j ( x ) β n ( x ′ ) , mathsf{Cov}(mathbf{x},mathbf{x}')approx hat{K}(mathbf{x},mathbf{x}')=sum_{j=1}^{N}sum_{n=1}^{N}k_{jn}alpha_{j}(mathbf{x})beta_{n}(mathbf{x}'), Cov(x,x′)≈K^(x,x′)=j=1∑Nn=1∑Nkjnαj(x)βn(x′),其中 k j n ∈ R k_{jn}inmathbb{R} kjn∈R为系数, α 1 ( x ) , … , α N ( x ) alpha_1(mathbf{x}),dots,alpha_{N}(mathbf{x}) α1(x),…,αN(x)和 β 1 ( x ′ ) , … , β N ( x ′ ) beta_1(mathbf{x}'),dots,beta_{N}(mathbf{x}') β1(x′),…,βN(x′)是线性无关的函数。 α j alpha_j αj和 β n beta_n βn的选择就有很多了。
Nyström methods
Nyström methods是把积分用数值的方式近似
∑
j
=
1
N
w
j
C
o
v
(
x
,
x
j
)
φ
^
i
(
x
j
)
=
λ
^
i
φ
^
i
(
x
)
.
sum_{j=1}^{N}w_jmathsf{Cov}(mathbf{x},mathbf{x}_j)hat{varphi}_i(mathbf{x}_j)=hat{lambda}_i hat{varphi}_i(mathbf{x}).
j=1∑NwjCov(x,xj)φ^i(xj)=λ^iφ^i(x). 进一步在积分点上解上式,即
∑
j
=
1
N
w
j
C
o
v
(
x
n
,
x
j
)
φ
^
i
(
x
j
)
=
λ
^
i
φ
^
i
(
x
n
)
,
n
=
1
,
…
,
N
.
sum_{j=1}^{N}w_jmathsf{Cov}(mathbf{x}_n,mathbf{x}_j)hat{varphi}_i(mathbf{x}_j)=hat{lambda}_ihat{varphi}_i(mathbf{x}_n),quad n=1,dots,N.
j=1∑NwjCov(xn,xj)φ^i(xj)=λ^iφ^i(xn),n=1,…,N.上式可以写成
C
W
y
i
=
λ
^
i
y
i
,
mathbf{C}mathbf{W}mathbf{y}_i=hat{lambda}_imathbf{y}_i,
CWyi=λ^iyi,
W
1
2
C
W
1
2
(
W
1
2
y
i
)
=
λ
^
i
W
1
2
y
i
,
mathbf{W}^{frac{1}{2}}mathbf{C}mathbf{W}^{frac{1}{2}}(mathbf{W}^{frac{1}{2}}mathbf{y}_i)=hat{lambda}_imathbf{W}^{frac{1}{2}}mathbf{y}_i,
W21CW21(W21yi)=λ^iW21yi,
C
mathbf{C}
C是对称半正定矩阵,
W
mathbf{W}
W是对角阵,
y
i
mathbf{y}_i
yi是向量。令
y
i
∗
=
W
1
2
y
i
mathbf{y}_i^{*}=mathbf{W}^{frac{1}{2}}mathbf{y}_i
yi∗=W21yi,解这个求特征向量和特征值的方程,注意特征向量的正规化。
φ
^
i
(
x
)
=
1
λ
^
i
∑
j
=
1
N
w
j
y
i
,
j
∗
C
o
v
(
x
,
x
j
)
.
hat{varphi}_i(mathbf{x})=frac{1}{hat{lambda}_i}sum_{j=1}^{N}sqrt{w_j}y_{i,j}^{*}mathsf{Cov}(mathbf{x},mathbf{x}_j).
φ^i(x)=λ^i1j=1∑Nwjyi,j∗Cov(x,xj).
Projection methods
残差 r IEVP ( x ) r_{text{IEVP}}(mathbf{x}) rIEVP(x): r IEVP ( x ) = ∑ j = 1 N d j i ( ∫ Ω C o v ( x , x ′ ) h j ( x ′ ) d x ′ − λ ^ i h j ( x ) ) . r_{text{IEVP}}(mathbf{x})=sum_{j=1}^{N}d_j^{i}left(int_{Omega}mathsf{Cov}(mathbf{x},mathbf{x}')h_j(mathbf{x}')mathsf{d}mathbf{x}'-hat{lambda}_i h_j(mathbf{x})right). rIEVP(x)=j=1∑Ndji(∫ΩCov(x,x′)hj(x′)dx′−λ^ihj(x)).
(1) Collocation methods
希望 r IEVP ( x ) = 0 r_{text{IEVP}}(mathbf{x})=0 rIEVP(x)=0给定一些点 { x l } l = 1 P {mathbf{x}_{l}}_{l=1}^{P} {xl}l=1P。构造 P × N Ptimes N P×N矩阵 A = ( a i j : = ∫ Ω C o v ( x i , x ) h j ( x ) d x ) mathbf{A}=(a_{ij}:=int_{Omega}mathsf{Cov}(mathbf{x}_i,mathbf{x})h_j(mathbf{x})mathsf{d}mathbf{x}) A=(aij:=∫ΩCov(xi,x)hj(x)dx), N = ( n i j : = h j ( x i ) ) mathbf{N}=(n_{ij}:=h_j(mathbf{x}_i)) N=(nij:=hj(xi))。构造向量 ( d i ) j = d j i (mathbf{d}_i)_j=d_j^{i} (di)j=dji A d i = λ ^ i N d i . mathbf{A}mathbf{d}_i=hat{lambda}_imathbf{N}mathbf{d}_i. Adi=λ^iNdi.解它。
(2) Galerki methods
希望残差能和
h
j
h_j
hj张成的线性子空间正交,所以希望解
∫
Ω
r
IEVP
(
x
)
h
j
(
x
)
d
x
=
0
,
∀
j
=
1
,
…
,
N
.
int_{Omega}r_{text{IEVP}}(mathbf{x})h_j(mathbf{x})mathsf{d}mathbf{x}=0,quad forall j=1,dots,N.
∫ΩrIEVP(x)hj(x)dx=0,∀j=1,…,N.
构造
N
×
N
Ntimes N
N×N矩阵
B
=
(
b
l
n
:
=
∫
Ω
h
l
(
x
)
∫
Ω
C
o
v
(
x
,
x
′
)
h
n
(
x
′
)
d
x
′
d
x
)
mathbf{B}=(b_{ln}:=int_{Omega}h_l(mathbf{x})int_{Omega}mathsf{Cov}(mathbf{x},mathbf{x}')h_n(mathbf{x}')mathsf{d}mathbf{x}'mathsf{d}mathbf{x})
B=(bln:=∫Ωhl(x)∫ΩCov(x,x′)hn(x′)dx′dx),
M
=
(
m
l
n
:
=
∫
Ω
h
l
(
x
)
h
n
(
x
)
d
x
)
mathbf{M}=(m_{ln}:=int_{Omega}h_{l}(mathbf{x})h_{n}(mathbf{x})mathsf{d}mathbf{x})
M=(mln:=∫Ωhl(x)hn(x)dx)。
B
d
i
=
λ
^
i
M
d
i
.
mathbf{B}mathbf{d}_i=hat{lambda}_imathbf{M}mathbf{d}_i.
Bdi=λ^iMdi.
论文中有更细致的关于解法的分析。这里就不继续探究了,主要是简单了解一下方法,以后用的时候知道应该往哪个方向查。
最后
以上就是酷酷跳跳糖为你收集整理的Karhunen-Loève expansion and random field的全部内容,希望文章能够帮你解决Karhunen-Loève expansion and random field所遇到的程序开发问题。
如果觉得靠谱客网站的内容还不错,欢迎将靠谱客网站推荐给程序员好友。
发表评论 取消回复