Learned IE: Solving the minimization problem¶
This is the first supplementary notebook for the paper [1] about learned infinite elements. It demonstrates how the minimization problem is solved and visualizes the resulting approximation of the DtN function. Two different kinds of exterior domains are treated:
- A) Homogeneous medium.
- B) Piecewise constant medium with a step discontinuity.
[1] Learned infinite elements (Hohage, Lehrenfeld, Preuß)
Most of the imports are standard except the module ceres_dtn, which provides the functionality for solving the minimization problem.
import numpy as np
from math import sqrt,exp
np.random.seed(123)
from scipy.special import hankel1,h1vp,jv,yv,jvp,yvp
import matplotlib.pyplot as plt
from matplotlib import gridspec
import ceres_dtn as opt
We set some graphic options to improve the presentation.
from IPython.core.display import display, HTML
display(HTML("<style>.container { width:80% !important; }</style>"))
plt.rc('legend', fontsize=16)
plt.rc('axes', titlesize=16)
plt.rc('axes', labelsize=16)
#%matplotlib notebook
%matplotlib inline
A) Homogeneous medium¶
In this problem a homogeneous medium in the exterior of a disk with radius $r = a$ is treated: $$ - \Delta u - k^2 u = 0, \quad r \geq a, $$ for some constant $k > 0$. The parameters are chosen below.
k = 16
a = 1.0
The main idea of our method is to learn an (efficient) exterior discretization which represents the DtN operator as good as possible. The reference DtN operator is assumed to diagonalize in the basis of the eigenfunctions of the Laplace-Beltrami operator on the coupling boundary. Also the eigenvalues of the continuous eigenpairs $$ - \Delta u_{\ell} = \lambda_{\ell} u_{\ell}, \quad \lambda_{\ell} = (\ell / a)^2 $$ are needed for implementation of the method. The diagonal elements of the DtN in this eigenbasis are given by $$ \text{dtn}^{\text{ref}}(\lambda) = - k \frac{ (H^{(1)}_{ a \sqrt{\lambda }})^{\prime}(k a) }{ (H^{(1)}_{a \sqrt{\lambda }})(k a) } , $$ for $\lambda = \lambda_{\ell}$ with the Hankel functions $H^{(1)}_{a \sqrt{\lambda } }$ of the first kind of order $ a \sqrt{\lambda}$. The continuous DtN function is defined below.
def dtn_ref(nu):
return -k*h1vp(nu,k*a) / hankel1(nu,k*a)
Next it has to be decided how many of the eigenvalues should be considered for approximating the continuous DtN function. This is controlled by the parameter $\ell = 0, \ldots, $ Lmax. Furthermore, each of the continuous DtN numbers is associated with a weight factor $ w_{\ell} \sim \exp(- \ell \alpha) $ that controls how much emphasis should be placed on fitting this particular mode. The decay of the weights is controlled by the parameter alpha_decay.
Lmax = 72
alpha_decay = 2/3
Now the arrays for the eigenvalues lam, the DtN numbers dtn_nr and the weights weights are set up. Additionally, we define an empty array final_res,
which will later be passed to the minimization routine to collect the final residuals.
lam = np.array([(l/a)**2 for l in range(Lmax)])
dtn_nr = np.array([ dtn_ref(sqrt(lami)*a) for lami in lam ])
weights = np.array([10**6*np.exp(-l*alpha_decay) for l in range(Lmax)])
final_res = np.zeros(len(lam),dtype=float)
These arrays are used to initialize an instance of the learned DtN class.
l_dtn = opt.learned_dtn(lam,dtn_nr,weights**2)
The instance l_dtN has a method called Run which solves the minimization problem
$$ \min_{A,B} \sum\limits_{\ell}{ \vert w_{\ell}( \text{dtn}^{\text{ref}}(\lambda_{\ell}) - \text{dtn}(\lambda_{\ell}) ) \vert^2}. $$
To this end, the method requires initial guesses A_guess, B_guess for the N+1x N+1 matrices $A$ and $B$, which are chosen as follows. For N = 0 we simply start with a random guess while for N > 1 the minimizer from step N-1 is recycled. This is implemented in the function new_initial_guessbelow.
def new_initial_guess(A_old,B_old,ansatz):
N = A_old.shape[0]
A_guess = np.zeros((N+1,N+1),dtype='complex')
B_guess = np.zeros((N+1,N+1),dtype='complex')
if ansatz in ["medium","full"]:
A_guess = 1e-3*(np.random.rand(N+1,N+1) + 1j*np.random.rand(N+1,N+1))
B_guess = 1e-3*(np.random.rand(N+1,N+1) + 1j*np.random.rand(N+1,N+1))
A_guess[:N,:N] = A_old[:]
B_guess[:N,:N] = B_old[:]
A_guess[N,N] = 1.0
elif ansatz == "minimalIC":
A_guess = 1e-3*(np.random.rand(N+1,N+1) + 1j*np.random.rand(N+1,N+1))
A_guess[:N,:N] = A_old[:]
B_guess[:N,:N] = B_old[:]
A_guess[N,N] = -100-100j
B_guess[0,N] = 1e-3*(np.random.rand(1) + 1j*np.random.rand(1))
return A_guess,B_guess
The minimization will be run for $N = 0, \ldots,$ Nmax-1. For later use we define some lists Lone,Ltwo and relative_residuals which are used to store the obtained matrices and relative residuals $ \vert \text{dtn}^{\text{ref}}(\lambda_{\ell}) - \text{dtn}(\lambda_{\ell}) \vert / \vert \text{dtn}^{\text{ref}}(\lambda_{\ell}) \vert $ for each $N$.
Furthermore, we have to define the ansatzfor $\text{dtn}$. The learned DtN function is in general of the form
$$ \text{dtn}(\lambda_{\ell}) = \left( S_{00}^{\ell} - \sum\limits_{i=1}^{N } \sum\limits_{j=1}^{ N } S_{0 i}^{\ell} \left[ \left( S^{\ell}_{EE} \right)^{-1} \right]_{i-1 ,j-1} S_{j 0}^{\ell} \right), \quad (*) $$
with $S^{\ell} := A +\lambda_{\ell} B$.
The following choices for ansatz are currently supported:
the complicated ansatz shown above for the minimization.
full: For this choice the matrices $A,B$ are full matrices with no specific structure. For this choice one has to work with
medium: In this case the block of the matrices that describes the exterior is assumed to be diagonal. This allows to work with the simplified ansatz $$ \text{dtn}(\lambda) = A_{00} + \lambda B_{00} - \sum\limits_{j=1}^{N}{ \frac{(A_{0j} + \lambda B_{0j}) ( A_{j0} + \lambda B_{j0} ) }{ A_{jj} + \lambda B_{jj} } }. $$minimalIC: A small additional reduction is obtained by setting some entries of $B$ equal to one $$ \text{dtn}(\lambda) = A_{00} + \lambda B_{00} - \sum\limits_{j=1}^{N}{ \frac{(A_{0j} + \lambda B_{0j}) ( A_{j0} + \lambda ) }{ A_{jj} + \lambda } }. $$
The behavior of the solver can be controlled by passing a dictionary called flags. The meaning of the options is described here.
Nmax = 7
A_IE = []
B_IE = []
relative_residuals = []
ansatz = "minimalIC"
#ansatz = "medium"
flags = {"max_num_iterations":5000000,
"use_nonmonotonic_steps":True,
"minimizer_progress_to_stdout":False,
"num_threads":4,
"report_level":"Brief",
"function_tolerance":1e-6,
"parameter_tolerance":1e-8}
A_guess = np.random.rand(1,1) + 1j*np.random.rand(1,1)
B_guess = np.random.rand(1,1) + 1j*np.random.rand(1,1)
Now the mininization loop over $N$ is run. During the call Run the initial guesses l1_guess,l2_guess are overwritten
with the obtained minimizers. For each $N$ the solver provides a short report about the number of required iterations as well as the
initial and final costs. The latter gives a good indicator to assess if the minimization has been successfull. If the solver should ever terminate with NO CONVERGENCE it is useful to set the option "report_level":"Full"
in the flags. A more detailed report is then provided which may help to identify the issue.
for N in range(Nmax):
l_dtn.Run(A_guess,B_guess,ansatz,flags,final_res)
A_IE.append(A_guess.copy()), B_IE.append(B_guess.copy()),relative_residuals.append(final_res.copy())
A_guess,B_guess = new_initial_guess(A_IE[N],B_IE[N],ansatz)
To check if the minimization has been successfull it makes sense to evaluate the learned dtn function dtn_learned
at some sample points lam_sample and compare to to the reference.
def dtn_learned(A_IE,B_IE,lami):
n_L,_ = A_IE.shape
S_lam = A_IE + lami*B_IE
result = S_lam[0,0]
if n_L > 1:
invS_lam = np.linalg.inv(S_lam[1:,1:])
for i in range(1,n_L):
for j in range(1,n_L):
result -= S_lam[0,i]*invS_lam[i-1,j-1]*S_lam[j,0]
return result
lam_sample = np.linspace(0.5,1000,num=200)
dtn_ref_sampled = np.array([ dtn_ref(sqrt(lami)*a) for lami in lam_sample])
Moreover, for a more qualitative check the relative residuals $ \vert \text{dtn}^{\text{ref}}(\lambda_{\ell}) - \text{dtn}(\lambda_{\ell}) \vert / \vert \text{dtn}^{\text{ref}}(\lambda_{\ell}) \vert, \ell = 0, \ldots, $ Lmax for different $N$ are also plotted in the left Figure. The real and imaginary part of the analytic DtN
are compared to the learned DtN function for different $N$ in the middle and the right Figure.
fig = plt.figure(figsize=(24, 8))
gs = gridspec.GridSpec(1, 3)
ax0 = plt.subplot(gs[0])
plt.title("Relative residuals")
plt.xlabel("$\ell$")
ax1 = plt.subplot(gs[1])
plt.xlabel("$\lambda$")
plt.title("Real part")
ax2 = plt.subplot(gs[2])
plt.xlabel("$\lambda$")
plt.title("Imaginary part")
ax1.plot(lam_sample, dtn_ref_sampled.real,linewidth=8,color='lightgray')
ax2.plot(lam_sample, dtn_ref_sampled.imag,label="ref",linewidth=8,color='lightgray')
for N in range(Nmax):
if N % 2 == 0:
dtn_approx = np.array([dtn_learned(A_IE[N],B_IE[N],lami) for lami in lam_sample])
ax0.semilogy(np.arange(Lmax),relative_residuals[N])
ax1.plot(lam_sample, dtn_approx.real)
ax2.plot(lam_sample, dtn_approx.imag,label="$N$ = {0}".format(N))
plt.legend()
The so callled poles are given by those complex numbers $\lambda$ for which $A_{EE}+\lambda B_{EE}$ is not invertible (here $E$ represents exterior degrees of freedom). Notice that at these $\lambda$ the formula (*) breaks down. Below these poles are computed as the negative eigenvalues of $ \left[ B_{EE} \right]^{-1} A_{EE}$ and plotted for the different $N$.
colors = ['b', 'c', 'y', 'm', 'r']
marker = ['*','o','s','x','.']
fig, ax = plt.subplots(figsize=(8,8))
for N,clr,mkr in zip(range(2,Nmax),colors,marker):
poles = -np.linalg.eigvals(np.linalg.inv(B_IE[N][1:,1:]) @ A_IE[N][1:,1:])
ax.scatter(poles.real,poles.imag,s=100, c=clr, label="$N$ = {0}".format(N),marker=mkr,
alpha=1.0, edgecolors='none')
ax.legend()
ax.grid(True)
plt.show()
B) Piecewise constant medium with a step discontinuity.¶
Let us get ready to tackle a more complicated problem.
from scipy.special import jv,yv,jvp,yvp
A_IE = []
B_IE = []
relative_residuals = []
This problem involves an inhomogeneous exterior domain. The wavenumber has a step discontinuity defined by:
$$ k(r) = \begin{cases} k_{I} & r < R_{\infty}, \\ k_{\infty} & r > R_{\infty}, \end{cases}$$ for some $ R_{\infty} > a$ where $a$ is again the radius of the coupling boundary. The parameters can be chosen below.
a = 1.0
R_inf = 2.0
k_I = 16
k_inf = 8
alpha_decay = 1.5
To determine the DtN numbers for each mode $\nu$ one has to solve
$$ -\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_{\nu} }{\partial r} \right) + \left( \omega^2(r) - \frac{\nu^2}{r^2} \right)u_{\nu} = 0 \quad r > a,$$ with boundary condition $u_{\nu}(a) = 1$ and radiation condition at infinity. The solution is of the form
$$ u_{\nu}(r) = \begin{cases} A_{\nu} J_{\nu}(k_{I}r) + B_{\nu} Y_{\nu}(k_I r) & r < R_{\infty}, \\ C_{\nu} H_{\nu}^{(1)}(k_{\infty}r) & r > R_{\infty}, \end{cases} $$
with constants $A_{\nu}, B_{\nu}$ and $C_{\nu}$ to be determined by the boundary condition at $r = a$ and continuity requirements at $r = R_{\infty}$. This leads to the system of equations
$$ A_{\nu} J_{\nu}(k_I a ) + B_{\nu} Y_{\nu}(k_I a ) = 1 .$$ $$ A_{\nu} J_{\nu}(k_I R ) + B_{\nu} Y_{\nu}(k_I R ) - C_{\nu} H_{\nu}^{(1)}(k_{\infty} R) = 0. $$ $$ A_{\nu} k_I J_{\nu}^{\prime}(k_I R ) + B_{\nu} k_I Y_{\nu}^{\prime}(k_I R ) - C_{\nu} k_{\infty} (H_{\nu}^{(1)})^{\prime}(k_{\infty} R) = 0. $$
A calculation shows that $$A_{\nu} = (1-B_{\nu} Y_{\nu}(k_I a )) / J_{\nu}(k_I a ), $$ $$ B_{\nu} = \frac{1}{M_{\nu}} \left[ \frac{k_{\infty}}{k_I} (H_{\nu}^{(1)})^{\prime}(k_{\infty} R) \frac{J_{\nu}(k_I R )}{ J_{\nu}(k_I a ) } - H_{\nu}^{(1)}(k_{\infty} R) \frac{ J_{\nu}^{\prime}(k_I R )}{ J_{\nu}(k_I a ) } \right], $$ for a certain $M_{\nu}$ defined in the code below.
def M_nu(nu):
tmp1 = -(k_inf/k_I)*(h1vp(nu,k_inf*R_inf)/jv(nu,k_I*a))*( jv(nu,k_I*a)*yv(nu,k_I*R_inf) - yv(nu,k_I*a)*jv(nu,k_I*R_inf) )
tmp2 = (hankel1(nu,k_inf*R_inf)/jv(nu,k_I*a))*( yvp(nu,k_I*R_inf)*jv(nu,k_I*a) - yv(nu,k_I*a)*jvp(nu,k_I*R_inf))
return tmp1+tmp2
def B_nu(nu):
tmp1 = (k_inf/k_I)*h1vp(nu,k_inf*R_inf)*jv(nu,k_I*R_inf)/jv(nu,k_I*a)
tmp2 = -hankel1(nu,k_inf*R_inf)*jvp(nu,k_I*R_inf)/jv(nu,k_I*a)
return (tmp1+tmp2)/M_nu(nu)
def A_nu(nu):
return (1-B_nu(nu)*yv(nu,k_I*a))/jv(nu,k_I*a)
def C_nu(nu):
tmp1 = ( yvp(nu,k_I*R_inf)*jv(nu,k_I*a) - yv(nu,k_I*a)*jvp(nu,k_I*R_inf) )*( jv(nu,k_I*R_inf) / jv(nu,k_I*a)**2)
tmp2 = -( jv(nu,k_I*a)*yv(nu,k_I*R_inf) - yv(nu,k_I*a)*jv(nu,k_I*R_inf) )*( jvp(nu,k_I*R_inf) / jv(nu,k_I*a)**2 )
return (tmp1+tmp2)/M_nu(nu)
We check that the coefficients defined above indeed solve the system
$$ \begin{bmatrix} J_{\nu}(k_I a ) & Y_{\nu}(k_I a ) & 0 \\ J_{\nu}(k_I R ) & Y_{\nu}(k_I R ) & - H_{\nu}^{(1)}(k_{\infty} R) \\ \omega_I J_{\nu}^{\prime}(k_I R ) & k_I Y_{\nu}^{\prime}(k_I R ) & - k_{\infty} (H_{\nu}^{(1)})^{\prime}(k_{\infty} R) \\ \end{bmatrix} \begin{bmatrix} A_{\nu} \\ B_{\nu} \\ C_{\nu} \end{bmatrix} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, $$ for $\nu = 0,1, \ldots$.
def Check(nu):
M = np.array( [ [jv(nu,k_I*a),yv(nu,k_I*a),0],
[jv(nu,k_I*R_inf),yv(nu,k_I*R_inf),-hankel1(nu,k_inf*R_inf)],
[k_I*jvp(nu,k_I*R_inf),k_I*yvp(nu,k_I*R_inf),-k_inf*h1vp(nu,k_inf*R_inf)]
])
print( "m = {0}, diff = {1}".format(nu, np.linalg.norm( M @ np.array([A_nu(nu),B_nu(nu),C_nu(nu)]) - np.array([1,0,0])) ))
for nu in range(20):
Check(nu)
The reference DtN function is then given by: $$ \text{dtn}^{\text{ref}}(\nu) = -\frac{\partial u_{\nu}(a)}{\partial r } = - k_I \left[ A_{\nu} J_{\nu}^{\prime}(k_I a ) + B_{\nu} Y_{\nu}^{\prime}(k_I a ) \right]. $$
def dtn_ref(nu):
return -k_I*( A_nu(nu)*jvp(nu,k_I*a) + B_nu(nu)*yvp(nu,k_I*a))
Similar as above we define the arrays lam, dtn_nr and weight, create an instance of the learned DtN and solve the
minimization problem using the Run method.
lam = np.array([(mm/a)**2 for mm in range(Lmax)])
dtn_nr = np.array([ -k_I*( A_nu(sqrt(lami)*a)*jvp(sqrt(lami)*a,k_I*a) + B_nu(sqrt(lami)*a)*yvp(sqrt(lami)*a,k_I*a)) for lami in lam ])
weights = np.array([10**6*exp(-mm/alpha_decay) for mm in range(len(lam))])
l_dtn = opt.learned_dtn(lam,dtn_nr,weights**2)
A_guess = np.random.rand(1,1) + 1j*np.random.rand(1,1)
B_guess = np.random.rand(1,1) + 1j*np.random.rand(1,1)
for N in range(Nmax):
l_dtn.Run(A_guess,B_guess,ansatz,flags,final_res)
A_IE.append(A_guess.copy()), B_IE.append(B_guess.copy()),relative_residuals.append(final_res.copy())
A_guess,B_guess = new_initial_guess(A_IE[N],B_IE[N],ansatz)
Let us again visualize the results.
lam_sample = np.linspace(0.5,800,num=10000)
dtn_ref_sampled = np.array([ dtn_ref(sqrt(lami)*a) for lami in lam_sample])
fig = plt.figure(figsize=(24, 8))
gs = gridspec.GridSpec(1, 3)
ax0 = plt.subplot(gs[0])
plt.title("Relative residuals")
plt.xlabel("$\ell$")
ax1 = plt.subplot(gs[1])
plt.xlabel("$\lambda$")
plt.title("Real part")
ax2 = plt.subplot(gs[2])
plt.xlabel("$\lambda$")
plt.title("Imaginary part")
ax1.plot(lam_sample, dtn_ref_sampled.real,linewidth=8,color='lightgray')
ax2.plot(lam_sample, dtn_ref_sampled.imag,label="ref",linewidth=8,color='lightgray')
for N in range(Nmax):
if N % 2 == 0:
dtn_approx = np.zeros(len(lam_sample),dtype=complex)
opt.eval_dtn_fct(A_IE[N],B_IE[N],lam_sample,dtn_approx)
ax0.semilogy(np.arange(Lmax),relative_residuals[N])
ax1.plot(lam_sample, dtn_approx.real)
ax2.plot(lam_sample, dtn_approx.imag,label="$N$ = {0}".format(N))
plt.legend(loc='upper right')
And also plot the corresponding poles.
colors = ['b', 'c', 'y', 'm', 'r']
marker = ['*','o','s','x','.']
fig, ax = plt.subplots(figsize=(8,8))
for N,clr,mkr in zip(range(3,Nmax),colors,marker):
poles = -np.linalg.eigvals(np.linalg.inv(B_IE[N][1:,1:]) @ A_IE[N][1:,1:])
ax.scatter(poles.real,poles.imag,s=100, c=clr, label="$N$ = {0}".format(N),marker=mkr,
alpha=1.0, edgecolors='none')
ax.legend()
ax.grid(True)
plt.show()
