Find the LCA (lowest common ancestor) of two nodes on a binary tree


Reproducing Kernel Hilbert Space

Since writing about the simulation of fractional Brownian motion, I’ll spend an hour to write RKHS.

Definition: A Hilbert space \mathcal{H} is a RKHS if |\mathcal{F}_t[f]|=|f(t)|\leq M\|f\|_{\mathcal{H}}, \forall f\in \mathcal{H}

Theorem: If \mathcal{H} is a RKHS, then for each t \in X there exists a function K_t \in \mathcal{H} (called the representer of t) with the reproducing property

\mathcal{F}_t[f]=<K_t,f>_{\mathcal{H}}=f(t), \forall f\in\mathcal{H}

Therefore, K_t(t^{'})=<K_t,K_{t^{'}}>_{\mathcal{H}}.

Definition: K: X\times X\rightarrow\mathbb{R} is a reproducing kernel if it’s symmetric and positive definite.

Theorem: A RKHS defines a corresponding reproducing kernel. Conversely, a reproducing kernel defines a unique RKHS.

Once we have the kernels, if f(\cdot)=\sum \alpha_i K(t_i,\cdot), g(\cdot)=\sum \beta_i K(t^{'}_i,\cdot),

then <f,g>_{\mathcal{H}}=\sum\sum \alpha_i \beta_j K(t_i, t^{'}_j)

When it comes to the fractional Brownian motion

Theorem: For fBm RKHS K(x,x^{'})=R(x,x^{'}) is symmetric and positive definite, , there exists K^H(x,\cdot) s.t K(x,x^{'})=\int K^H(x,y)K^H(x^{'},y)dy

Take Wiener integral as an example:

\begin{pmatrix}    \textit{dual space}&\textit{inner product}&E&\leftrightarrow&<\cdot,\cdot>_{\mathcal{H}}&\leftrightarrow&<\cdot,\cdot>_{L^{2}}\\    \delta_{t}(\cdot)& &B^H_t&\leftrightarrow&R(\cdot,t)&\leftrightarrow&K^H(t,\cdot)\\    f(\cdot)& &\textit{"}\int_0^T f(t)B^H_tdt\textit{"}&\leftrightarrow&K^H\circ(K^H)^{*}f&\leftrightarrow&(K^H)^{*}f    \end{pmatrix}

where K^H\circ f(t)=\int_0^T K^H(t,s)f(s)ds,

(K^H)^{*}\circ f(t)=\int_0^T K^H(s,t)f(s)ds and

K^H\circ(K^H)^{*}f(t)=\int_0^T R(t,s)f(s)ds

For example, when H=\frac{1}{2}, K^H(t,s)=1_{[0,t]}(s), , we can check

E(B_tB_s)=<t\wedge\cdot, s\wedge\cdot>_{\mathcal{H}}=<1_{[0,t]}(\cdot),1_{[0,s]}(\cdot)>_{L^2}

Claim: By previous theorem, the reproducing kernel R(t,s) uniquely defines a RKHS L(R(t,\cdot)) with the inner product \mathcal{H}

how to simulate a fractional Brownian motion?

A fractional Brownian motion B^H_tis a Gaussian process satisfying

a) B^H_0=0

b) B^H_t is Gaussian

c) E(B^H_tB^H_s)=R(t,s)=\frac{1}{2}(t^{2H}+s^{2H}-|t-s|^{2H})

As for simulation, we just need to simulate the correlated Gaussian (\Delta B^H_{t_1},\dots,\Delta B^H_{t_n})\sim N(0,\Sigma), where

t_i=i*T/n=i \Delta t,

\Delta B^H_{t_i}=B^H_{t_i}-B^H_{t_{i-1}}, i=1,\dots,n,

\Sigma_{i,j}=E(B^H_{t_i}B^H_{t_j})=\frac{\Delta t^{2H}}{2}((i-j+1)^{2H}+(i-j-1)^{2H}-2(i-j)^{2H})

Now the question is just simulating L'\cdot N(0,I_{n\times n}) where L is the Cholesky decomposition of \Sigma, i.e (\Sigma=L'L)

Method 1: find \Sigma first and then use Cholesky decomposition to compute L. This works theoretically, but not practically because it’s too time and space costy to record \Sigma, not to mention decomposition step.

Method 2:

We denote \lambda(k)=\frac{(\Delta t)^{2H}}{2}((k+1)^{2H}+(k-1)^{2H}-2k^{2H})=\Sigma_{i,i+k}

Denote a circulant matrix

C=\begin{pmatrix}\lambda(0)&\lambda(1)&\lambda(2)&\dots&\lambda(n-1)&\lambda(n-2)&\dots&\lambda(1)\\ \lambda(1)&\lambda(0)&\lambda(1)&\dots&\lambda(n-2)&\lambda(n-1)&\dots&\lambda(2)\\ \lambda(2)&\lambda(1)&\lambda(0)&\dots&\lambda(n-3)&\lambda(n-2)&\dots&\lambda(3)\\ \vdots&\vdots&\vdots&\ddots&\dots&\dots&\dots&\dots\\    \lambda(n-1)&\lambda(n-2)&\lambda(n-3)&\dots&\lambda(0)&\lambda(1)&\dots&\lambda(n-2)\\    \dots    \end{pmatrix}_{2n-2\times 2n-2}

of which \Sigma is the left up n\times n block matrix.

The idea is, suppose C=DD^*=(D_1+iD_2)*(D_1^T-i D_2^T)where D_1, D_2 \in \mathbb{R}^{2n-2\times 2n-2}, now if we generate i.i.d \quad \epsilon^1,\epsilon^2\sim N(0,I_{2n-2,2n-2}),

We can conclude that

denote X_1+iX_2=D*(\epsilon^1+i\epsilon^2)=D_1\epsilon^1-D_2\epsilon^2+i(D_2\epsilon^1+D_1\epsilon^2)\in \mathbb{R}^{2n-2,1}, then X_1[1:n]\sim N(0,\Sigma) and X_2[1:n]\sim N(0,\Sigma)

let Q is the unitary matrix defined by

Q_{j,k}=(2n-2)^{-1/2}exp(-2i\pi\frac{jk}{2n-2}) and thus D=Q*\Lambda^{\frac{1}{2}}

so how to achieve \Lambda? The property of circulant matrix shows that

\lambda_k=\sum_{j=1}^{2n-2}c_j exp(-2i\frac{jk}{2n-2}), i.e

\Lambda=Q*\begin{pmatrix}\lambda(0)&\lambda(1)&\lambda(2)&\dots&\lambda(n-1)&\lambda(n-2)&\dots&\lambda(1)\\    \end{pmatrix}^{'} that can be achieved by FFT.

Next, according to the idea stated before, FFT(diag(\Lambda^{\frac{1}{2}})*(\epsilon^1+i\epsilon^2))=Q*diag(\Lambda^{\frac{1}{2}})*(\epsilon^1+i\epsilon^2)=X_1+iX_2

And the first n component of X_1 will be N(0,\Sigma)

Reference : Coeurjolly, JF. Simulation and identification of the fractional Brownian motion

Matlab code:

function []=fBM(n,H)
r=nan(n+1,1); r(1) = 1;
for k=1:n
r(k+1) = 0.5*((k+1)^(2*H) – 2*k^(2*H) + (k-1)^(2*H));
r=[r; r(end-1:-1:2)]; % first rwo of circulant matrix
lambda=real(fft(r))/(2*n); % eigenvalues
W = n^(-H)*cumsum(real(W(1:n+1))); % rescale
%Terminal=W(n+1);%terminal value of fBM