Spectrum of the preconditioned Hessian misfit operator

The linear source inversion problem

We consider the following linear source inversion problem. Find the state $u \in \mathcal{V} := H^1_{\Gamma_D}(\Omega)$ and the source (parameter) $m \in \mathcal{M} := H^1(\Omega)$ that solves

$$ \begin{aligned} {} & \min_m \underbrace{\frac{1}{2} \| Bu - \boldsymbol{d} \|^2 + \frac{1}{2} \int_\Omega \left[ \delta|m-m_0|^2 + \gamma|\nabla (m - m_0)|^2 \right] dx}_{J(u(m), m)} & {}\\ {\rm s.t.} & {} &{} \\ {} & -{\rm div}(k \nabla u) + cu = m & {\rm in} \; \Omega\\ {} & u = 0 & {\rm on } \; \Gamma_D\\ {} & k \frac{\partial u}{\partial n} = 0 & {\rm on } \; \Gamma_N\\ \end{aligned} $$

Here:

  • $\boldsymbol{d}$ is a $n_{\rm obs}$ finite dimensional vector that denotes noisy observations of the state $u$ in $n_{\rm obs}$ locations $\mathbf{x}_i$, $i=1,\ldots,n_{\rm obs}$. Specifically, $d_i = u_{\rm true}( {\bf x}_i ) + \eta_i$, where $\eta_i$ are i.i.d. $\mathcal{N}(0, \sigma^2)$.

  • $B: \mathcal{V} \rightarrow \mathbb{R}^{n_{\rm obs}}$ is the linear operator that evaluates the state $u$ at the observation locations $\mathbf{x}_i$, $i=1,\ldots,n_{\rm obs}$.

  • $\| \cdot \|$ denotes the Euclidean norm in $\mathbb{R}^{n_{\rm obs}}$.

  • $\delta$ and $\gamma$ are the parameters of the regularization penalizing the $L^2(\Omega)$ and $H^1(\Omega)$ norm of $m-m_0$, respectively.

  • $k$ and $c$ are given coefficients representing the diffusivity coefficient and the reaction term, respectively.

  • $\Gamma_D \subset \partial \Omega$, $\Gamma_N \subset \partial \Omega$ represents the subdomain of $\partial\Omega$ where we impose Dirichlet or Neumann boundary conditions, respectively.

The variational (or weak) form of the forward problem:

Find $u\in \mathcal{V}$ such that

$$ \int_{\Omega} k \nabla u \cdot \nabla \tilde{p} \, dx +\int_{\Omega} c\, u\,\tilde{p} \, dx = \int_{\Omega} m\,\tilde{p} \, dx, \text{ for all } \tilde{p} \in \mathcal{V}.$$

Gradient evaluation:

The Lagrangian functional $\mathscr{L}: \mathcal{V} \times \mathcal{M} \times \mathcal{V} \rightarrow \mathbb{R}$ is given by

$$ \mathscr{L}(u,m,p):= \underbrace{\frac{1}{2} \| Bu - \boldsymbol{d} \|^2 + \frac{1}{2} \int_\Omega \left[ \delta|m-m_0|^2 + \gamma|\nabla (m - m_0)|^2 \right] dx}_{J(u,m)} + \underbrace{\int_{\Omega} k \nabla u \cdot \nabla p \, dx +\int_{\Omega} c\, u\,p \, dx - \int_{\Omega} m\,p \, dx}_{\text{forward problem}}. $$

Then the gradient of the cost functional $\mathcal{J}(u(m), m)$ with respect to the parameter $m$ is

$$ (\mathcal{G}(m), \tilde m) := \mathscr{L}_m(u,m,p)(\tilde{m}) = \gamma \int_\Omega \nabla (m-m_0) \cdot \nabla \tilde{m}\, dx + \delta \int_{\Omega} (m-m_0) \cdot \tilde{m} \, dx - \int_\Omega \tilde{m}\, p\, dx \quad \forall \tilde{m} \in \mathcal{M}, $$

where $u \in \mathcal{V}$ is the solution of the forward problem,

$$ (\mathscr{L}_p(u,m,p), \tilde{p}) := \int_{\Omega} k \nabla u \cdot \nabla \tilde{p} \, dx +\int_{\Omega} c\, u\,\tilde{p} \, dx - \int_{\Omega} m\,\tilde{p} \, dx = 0 \quad \forall \tilde{p} \in \mathcal{V}, $$

and $p \in \mathcal{V}$ is the solution of the adjoint problem,

$$ (\mathscr{L}_u(u,m,p), \tilde{u}) := \int_{\Omega} k \nabla \tilde{u} \cdot \nabla p \, dx +\int_{\Omega} c\, \tilde{u} \,p \, dx + \int_\Omega \tilde{u}B^*(Bu - \boldsymbol{d})\,dx = 0 \quad \forall \tilde{u} \in \mathcal{V}.$$

Above $B^*: \mathbb{R}^{n_{\rm obs}} \mapsto \mathcal{V}$ is the adjoint of $B$, i.e. the operator defined as

$$ \underbrace{y^T (B u)}_{\text{inner product in } \mathbb{R}^{n_{\rm obs}}} = \underbrace{\int_\Omega u\, B^*y \, dx}_{\text{inner product in }\mathcal{V}} \quad forall y \in \mathbb{R}^{n_{\rm obs}},u \in \mathcal{V}. $$

Hessian action:

To evaluate the action $\mathcal{H}(m)(\hat{m})$ of the Hessian is a given direction $\hat{m}$, we consider variations of the meta-Lagrangian functional

$$ \begin{aligned} \mathscr{L}^H(u,m,p; \hat{u}, \hat{m}, \hat{p}) := & {} & {} \\ {} & \gamma \int_\Omega \nabla (m-m_0) \cdot \nabla \hat{m}\, dx + \delta \int_{\Omega} (m-m_0) \cdot \hat{m} \, dx - \int_\Omega \hat{m}\, p\, dx & \text{gradient}\\ {} & + \int_{\Omega} k \nabla u \cdot \nabla \hat{p} \, dx +\int_{\Omega} c\, u\,\hat{p} \, dx - \int_{\Omega} m\,\hat{p} \, dx & \text{forward eq}\\ {} & + \int_{\Omega} k \nabla \hat{u} \cdot \nabla p \, dx +\int_{\Omega} c\, \hat{u} \,p \, dx + \int_\Omega \hat{u}\,B^*(Bu - \boldsymbol{d})\,dx & \text{adjoint eq}. \end{aligned} $$

Then action of the Hessian is a given direction $\hat{m}$ is

$$ \begin{aligned} (\tilde{m}, \mathcal{H}(m)(\hat{m}) ) & := \mathscr{L}^H_m(u,m,p; \hat{u}, \hat{m}, \hat{p})(\tilde{m}) \\ {} & = \gamma \int_\Omega \nabla \tilde{m}\cdot \nabla \hat{m}\, dx + \delta \int_{\Omega} \tilde{m}\cdot \hat{m} \, dx - \int_{\Omega} \tilde{m} \,\hat{p} \, dx \quad \forall \tilde{m} \in \mathcal{M}, \end{aligned} $$

where

  • $u\in \mathcal{V}$ and $p \in \mathcal{V}$ are the solution of the forward and adjoint problem, respectively;

  • $\hat{u} \in \mathcal{V}$ is the solution of the incremental forward problem,

$$ \left( \mathscr{L}^H_p(u,m,p; \hat{u}, \hat{m}, \hat{p}),\tilde{p}\right) := \int_{\Omega} k \nabla \hat{u} \cdot \nabla \tilde{p} \, dx +\int_{\Omega} c\, \hat{u}\,\tilde{p} \, dx - \int_{\Omega} \hat{m}\,\tilde{p} \, dx = 0 \quad \forall \tilde{p} \in \mathcal{V}; $$
  • and $\hat{p} \in \mathcal{V}$ is the solution of the incremental adjoint problem, $$ \left( \mathscr{L}^H_u(u,m,p; \hat{u}, \hat{m}, \hat{p}), \tilde{u}\right) := \int_{\Omega} k \nabla \tilde{u} \cdot \nabla \hat{p} \, dx +\int_{\Omega} c\, \tilde{u} \,\hat{p} \, dx + \int_\Omega \tilde{u}B^*B\hat{u}\,dx = 0 \quad \forall \tilde{u} \in \mathcal{V}. $$

It worth to notice that the Hessian expression does not depend on the data $u_d$ or on $u$, $m$, $p$. This is a characteristic of all linear inverse problems.

Solution of inverse problem

Note that due to linearity of the inverse problem, Newton's method will converge in a single iteration. That is the solution of the inverse problem $m^*$ is the solution of the linear system $$ \left(\tilde{m}, \mathcal{H}\,m^* \right) = -\left(\tilde{m}, \mathcal{G}(0)\right) \quad \forall \tilde{m} \in \mathcal{M}, $$

where the evaluation of the gradient $\mathcal{G}(0)$ involve the solution $u$ and $p$ of the forward and adjoint problem (respectively) for $m = 0$. Similarly, the Hessian action in an arbitrary direction $\hat{m}$ requires to additional solve the incremental forward and adjoint problems.

Discrete optimality conditions:

$ \def\tu{\tilde u} \def\tm{\tilde m} \def\tp{\tilde p} \def\hu{\hat u} \def\hp{\hat p} \def\hm{\hat m} $ $ \def\bu{{\bf u}} \def\bm{{\bf m}} \def\bp{{\bf p}} \def\btu{{\bf \tilde u}} \def\btm{{\bf \tilde m}} \def\btp{{\bf \tilde p}} \def\bhu{{\bf \hat u}} \def\bhm{{\bf \hat m}} \def\bhp{{\bf \hat p}} \def\bg{{\bf g}} $ $ \def\bA{{\bf A}} \def\bB{{\bf B}} \def\bC{{\bf C}} \def\bH{{\bf H}} \def\bR{{\bf R}} \def\bW{{\bf W}} $

Let us introduce the finite dimensional subspace $\mathcal{M}_h \subset \mathcal{M}$ and $\mathcal{V}_h \subset \mathcal{V}$. Let $\{ \phi_i(\boldsymbol{x}) \}_{i=1}^{n_M}$ be a basis of $\mathcal{M}_h$, and $\{ \psi_i(\boldsymbol{x}) \}_{i=1}^{n_V}$ be a basis of $\mathcal{V}_h$.

Define the following matrices:

  • $\bA \in \mathbb{R}^{n_V \times n_V}$ stemming from Galerkin discretization of the left hand side of the forward problem, i.e. whose (i,j)-entries is given by $$ A_{ij} = \int_{\Omega} k \nabla \psi_j \cdot \nabla \psi_i \, dx +\int_{\Omega} c\, \psi_j \,\psi_i \, dx \quad \forall i,j=1,\ldots, n_V;$$

  • $\bC \in \mathbb{R}^{n_V \times n_M}$ stemming from Galerkin discretization of the right hand side of the forward problem, i.e. whose (i,j)-entries is given by

$$ C_{ij} = \int_{\Omega} \phi_j \,\psi_i \, dx \quad \forall i=1,\ldots, n_V, \forall j=1,\ldots, n_M$$

  • $\bB \in \mathbb{R}^{n_{\rm obs} \times n_V}$ stemming from the discretization of the observation operator $B$

  • $\bR$ stemming from Galerkin discretization of the Hessian of the regularization term, i.e. whose (i,j)-entries is given by $$ R_{ij} = \gamma \int_\Omega \nabla \phi_j \cdot \nabla \phi_i \, dx + \delta \int_{\Omega} \phi_j \, \phi_i \, dx $$

Let us also denote the vectors corresponding to the discretization of the functions $u, m, p$ by $\bu, \bm, \bp$.

Then, using the above notation, we can rewrite the inverse problem as

$$ \begin{aligned} {} & \min_\bm \frac{1}{2} \| \bB\,\bu - \boldsymbol{d} \|^2 + \frac{1}{2} (\bm-\bm_0)^T \bR (\bm-\bm_0) & {}\\ {\rm s.t.} & {} &{} \\ {} & \bA\,\bu = \bC \bm, \end{aligned} $$

where $\bm_0$ is the vector of degree of freedoms of the interpolant of $m_0$ onto $\mathcal{M}_h$.

The optimalitity conditions for the above discrete system can be readily obtained by eliminating the state $\bu = \bA^{-1}\bC \bm$.

In fact, first order optimality conditions for the mimization problem
$$ \min_\bm \frac{1}{2} \| \bB\bA^{-1}\bC \bm - \boldsymbol{d} \|^2 + \frac{1}{2} (\bm-\bm_0)^T \bR (\bm-\bm_0) $$ are $$ \underbrace{(\bC^T \bA^{-T}\bB^T\bB\bA^{-1}\bC +\bR)}_{\text{Hessian}} \bm = \underbrace{\bC^T \bA^{-T}\bB^T\boldsymbol{d} + \bR \bm_0}_{\text{gradient at }\bm = \boldsymbol{0} }.$$

Note that: 1) Computing the gradient at $\bm = \boldsymbol{0}$ requires solving

  • the forward problem $\bA\,\bu = \bC \boldsymbol{0}$, which has the trivial solution $\bu = \boldsymbol{0}$;
  • the adjoint problem $\bA^T \bp = -\bB^T\boldsymbol{d}$. 2) Each Hessian-vector multiplication in a direction $\hat{\bm}$ requires solving
  • One incremental forward problem $\bA \hat{\bu} = \bC \hat{\bm}$;
  • One incremental adjoint problem $\bA^T\hat{\bp} = -\bB^T\bB \hat{\bu}$.

Spectrum of the Hessian

We decompose the Hessian $\bH$ in two components

$$ \bH = \underbrace{\bC^T \bA^{-T}\bB^T\bB\bA^{-1}\bC}_{\text{Hessian of the data misfit} \bH_m} + \underbrace{\bR}_{\text{Hessian of the regularization}}. $$

We then consider the generalized eigenproblem

$$ \bH_m \hat{\bm}_i = \lambda_i \bR \hat{\bm}_i$$

Note that: 1) Eigenvectors $\hat{\bm}_i$ corresponding to large eigenvalues ($\lambda_i \gg 1$) denotes modes (directions) in parameter space that are strongly informed by the data and are weakly penalized by the regularization. 2) Eigenvectors $\hat{\bm}_i$ corresponding to small eigenvalues ($\lambda_i \ll 1$) denotes modes (directions) in parameter space that are not informed by the data and are strongly penalized by the regularization.

For this reason, we can refer to the number of eigenvalues larger than one as the information dimension which is always smaller than both the parameter and the data dimensions.

Scalable algorithms for the solution of inverse problems should have a cost---measured in terms of forward/adjoint/incremental solves---that is independent of parameter/data dimension and only depends on the information dimension.

1. Load modules

In [1]:
import matplotlib.pyplot as plt
%matplotlib inline

import dolfin as dl
import ufl
import numpy as np

import hippylib as hp
from hippylib import nb


import logging
logging.getLogger('FFC').setLevel(logging.WARNING)
logging.getLogger('UFL').setLevel(logging.WARNING)
dl.set_log_active(False)

2. The linear source inversion problem

In [2]:
def A_varf(u,p):
    return k*ufl.inner(ufl.grad(u), ufl.grad(p))*ufl.dx \
           + c*u*p*ufl.dx
        
def f_varf(p):
    return dl.Constant(0.)*p*ufl.dx
        
def C_varf(m,p):
           return -m*p*ufl.dx
    
def R_varf(m_trial, m_test, gamma, delta):
    return dl.Constant(delta)*m_trial*m_test*ufl.dx \
           + dl.Constant(gamma)*ufl.inner(ufl.grad(m_trial), ufl.grad(m_test))*ufl.dx

def u_boundary(x, on_boundary):
    return on_boundary and x[1] < dl.DOLFIN_EPS

class Hessian:
    def __init__(self, A_solver, At_solver, C,B,R):
        self.misfit_only = False
        self.A_solver = A_solver
        self.At_solver = At_solver
        self.C = C
        self.B = B
        self.R = R
        
        self.uhat = dl.Vector()
        self.C.init_vector(self.uhat,0)
        self.fwd_rhs = dl.Vector()
        self.C.init_vector(self.fwd_rhs,0)
        self.adj_rhs = dl.Vector()
        self.C.init_vector(self.adj_rhs,0)
        self.phat = dl.Vector()
        self.C.init_vector(self.phat,0)
        
    def init_vector(self, x, dim):
        self.R.init_vector(x,dim)
        
    def mult(self, x, y):
        y.zero()
        self.C.mult(x, self.fwd_rhs)
        self.A_solver.solve(self.uhat, self.fwd_rhs)
        self.B.transpmult(self.B*self.uhat, self.adj_rhs)
        self.At_solver.solve(self.phat, self.adj_rhs)
        self.C.transpmult(self.phat, y)
        
        if not self.misfit_only:
            y.axpy(1., self.R*x)
        

def solve(nx, ny, targets, gamma, delta, verbose=True):
    
    rand_gen = hp.Random()
    mesh = dl.UnitSquareMesh(nx, ny)
    Vh1 = dl.FunctionSpace(mesh, 'Lagrange', 1)
    
    Vh = [Vh1, Vh1, Vh1]
    if verbose:
        print("Number of dofs: STATE={0}, PARAMETER={1}, ADJOINT={2}".format(
            Vh[hp.STATE].dim(), Vh[hp.PARAMETER].dim(), Vh[hp.ADJOINT].dim()) )


    u_bdr = dl.Constant(0.0)
    bc = dl.DirichletBC(Vh[hp.STATE], u_bdr, u_boundary)
    
    u_trial = dl.TrialFunction(Vh[hp.STATE])
    u_test  = dl.TestFunction(Vh[hp.STATE])
    m_trial = dl.TrialFunction(Vh[hp.PARAMETER])
    m_test  = dl.TestFunction(Vh[hp.PARAMETER])
    
    A,f = dl.assemble_system(A_varf(u_trial, u_test), f_varf(u_test), bcs=bc )
    A_solver = hp.PETScLUSolver(mesh.mpi_comm())
    A_solver.set_operator(A)
    
    At,dummy = dl.assemble_system(dl.adjoint( A_varf(u_trial, u_test) ), dl.Constant(0.)*u_test*dl.dx, bcs=bc )
    At_solver = hp.PETScLUSolver(mesh.mpi_comm())
    At_solver.set_operator(At)
    
    C = dl.assemble(C_varf(m_trial, u_test))
    bc.zero(C)
    
    R = dl.assemble(R_varf(m_trial, m_test, gamma, delta))
    R_solver = dl.PETScLUSolver()
    R_solver.set_operator(R)

    mtrue_str = 'min(0.5,exp(-100*(pow(x[0]-0.35,2) +  pow(x[1]-0.7,2))))'
    mtrue = dl.interpolate(
        dl.Expression(mtrue_str, degree=5), Vh[hp.PARAMETER]).vector()
    m0 = dl.interpolate(dl.Constant(0.0), Vh[hp.PARAMETER]).vector()
     
    if verbose:
        print( "Number of observation points: {0}".format(targets.shape[0]) )
    
    B = hp.assemblePointwiseObservation(Vh[hp.STATE], targets)
               
    #Generate synthetic observations
    utrue = dl.Function(Vh[hp.STATE]).vector()
    A_solver.solve(utrue, -(C*mtrue) )
    d = B*utrue
    MAX = d.norm("linf")
    rel_noise = 0.01
    noise_std_dev = rel_noise * MAX
    rand_gen.normal_perturb(noise_std_dev, d)

    u = dl.Function(Vh[hp.STATE]).vector()
    m = m0.copy()
    p = dl.Function(Vh[hp.ADJOINT]).vector()
    mg = dl.Function(Vh[hp.PARAMETER]).vector()
    rhs_adj = dl.Function(Vh[hp.STATE]).vector()
    
    # Solve forward:
    A_solver.solve(u, -(C*m) )
    # rhs for adjoint
    B.transpmult(d-(B*u), rhs_adj)
    # solve adj problem
    At_solver.solve(p, rhs_adj)
    #gradient
    C.transpmult(p, mg)
    mg.axpy(1., R*m)
   

    H = Hessian(A_solver,At_solver,C,B,R)

    solver = hp.CGSolverSteihaug()
    solver.set_operator(H)
    solver.set_preconditioner( R_solver )
    solver.parameters["print_level"] = -1
    solver.parameters["rel_tolerance"] = 1e-9
    solver.solve(m, -mg)

    if solver.converged:
        if verbose:
            print ("CG converged in ", solver.iter, " iterations.")
    else:
        print( "CG did not converged.")
        raise

    # Solve forward:
    A_solver.solve(u, -(C*m) )
 
    if verbose:
        plt.figure(figsize=(18,8))
        plt.subplot(2, 3, 1)
        dl.plot(dl.Function(Vh[hp.PARAMETER], mtrue), title = "True source")
        plt.subplot(2, 3, 2)
        dl.plot(dl.Function(Vh[hp.STATE], utrue), title="True state")
        plt.subplot(2, 3, 3)
        nb.plot_pts(targets, d,mytitle="Observations")
        plt.subplot(2, 3, 4)
        dl.plot(dl.Function(Vh[hp.PARAMETER], m), title = "Reconstructed source")
        plt.subplot(2, 3, 5)
        dl.plot(dl.Function(Vh[hp.STATE], u), title="Reconstructed state")
        plt.subplot(2, 3, 6)
        nb.plot_pts(targets, B*u-d,mytitle="Misfit")
        plt.show()

    H.misfit_only = True
    k_evec = 80
    p_evec = 5
    if verbose:
        print ("Double Pass Algorithm. Requested eigenvectors: {0}; Oversampling {1}.".format(k_evec,p_evec))
    
    Omega = hp.MultiVector(m, k_evec+p_evec)
    rand_gen.normal(1., Omega)
    lmbda, U = hp.doublePassG(H, R, R_solver, Omega, k_evec)

    if verbose:
        plt.figure()
        nb.plot_eigenvalues(lmbda, mytitle="Generalized Eigenvalues")
        nb.plot_eigenvectors(Vh[hp.PARAMETER], U, mytitle="Eigenvectors", which=[0,1,2,5,10,15])
        plt.show()
        
    return lmbda, U, Vh[hp.PARAMETER], solver.iter

3. Solution of the source inversion problem

In [3]:
ndim = 2
nx = 32
ny = 32

ntargets = 256
np.random.seed(seed=1)
targets = np.random.uniform(0.1,0.9, [ntargets, ndim] )


gamma = 1e-5
delta = 1e-9

k = dl.Constant(1.0)
c = dl.Constant(0.1)

lmbda, U, Vm, nit = solve(nx,ny, targets, gamma, delta)

Number of dofs: STATE=1089, PARAMETER=1089, ADJOINT=1089
Number of observation points: 256
CG converged in  38  iterations.

Double Pass Algorithm. Requested eigenvectors: 80; Oversampling 5.

4. Mesh independence of the spectrum of the preconditioned Hessian misfit

In [4]:
gamma = 1e-5
delta = 1e-9

k = dl.Constant(1.0)
c = dl.Constant(0.1)

n = [16,32,64]
lmbda1, U1, Vm1, niter1 = solve(n[0],n[0], targets, gamma, delta,verbose=False)
lmbda2, U2, Vm2, niter2 = solve(n[1],n[1], targets, gamma, delta,verbose=False)
lmbda3, U3, Vm3, niter3 = solve(n[2],n[2], targets, gamma, delta,verbose=False)

print ("Number of Iterations: ", niter1, niter2, niter3)
plt.figure(figsize=(18,4))
nb.plot_eigenvalues(lmbda1, mytitle="Eigenvalues Mesh {0} by {1}".format(n[0],n[0]), subplot_loc=131)
plt.ylim([1e-3,1e11])
nb.plot_eigenvalues(lmbda2, mytitle="Eigenvalues Mesh {0} by {1}".format(n[1],n[1]), subplot_loc=132)
plt.ylim([1e-3,1e11])
nb.plot_eigenvalues(lmbda3, mytitle="Eigenvalues Mesh {0} by {1}".format(n[2],n[2]), subplot_loc=133)
plt.ylim([1e-3,1e11])

nb.plot_eigenvectors(Vm1, U1, mytitle="Mesh {0} by {1} Eigenvector".format(n[0],n[0]), which=[0,1,5])
nb.plot_eigenvectors(Vm2, U2, mytitle="Mesh {0} by {1} Eigenvector".format(n[1],n[1]), which=[0,1,5])
nb.plot_eigenvectors(Vm3, U3, mytitle="Mesh {0} by {1} Eigenvector".format(n[2],n[2]), which=[0,1,5])

plt.show()

Number of Iterations:  37 38 37

5. Dependence on regularization parameters

We solve the problem for different values of the regularization parameters.

In [5]:
gamma = [1e-4, 1e-5, 1e-6]
delta = [1e-8, 1e-9, 1e-10]

k = dl.Constant(1.0)
c = dl.Constant(0.1)

lmbda1, U1, Vm1, niter1 = solve(nx,ny, targets, gamma[0], delta[0],verbose=False)
lmbda2, U2, Vm2, niter2 = solve(nx,ny, targets, gamma[1], delta[1],verbose=False)
lmbda3, U3, Vm3, niter3 = solve(nx,ny, targets, gamma[2], delta[2],verbose=False)

print ("Number of Iterations: ", niter1, niter2, niter3)
plt.figure(figsize=(18,4))
nb.plot_eigenvalues(lmbda1, mytitle="Eigenvalues gamma={0:1.1e}".format(gamma[0]), subplot_loc=131)
plt.ylim([1e-3,1e12])
nb.plot_eigenvalues(lmbda2, mytitle="Eigenvalues gamma={0:1.1e}".format(gamma[1]), subplot_loc=132)
plt.ylim([1e-3,1e12])
nb.plot_eigenvalues(lmbda3, mytitle="Eigenvalues gamma={0:1.1e}".format(gamma[2]), subplot_loc=133)
plt.ylim([1e-3,1e12])

nb.plot_eigenvectors(Vm1, U1, mytitle="gamma={0:1.1e} Eigenvector".format(gamma[0]), which=[0,1,5])
nb.plot_eigenvectors(Vm2, U2, mytitle="gamma={0:1.1e} Eigenvector".format(gamma[1]), which=[0,1,5])
nb.plot_eigenvectors(Vm3, U3, mytitle="gamma={0:1.1e} Eigenvector".format(gamma[2]), which=[0,1,5])

plt.show()

Number of Iterations:  23 38 64

6. Dependence on the PDE coefficients

Assume a constant reaction term $c = 0.1$, and we consider different values for the diffusivity coefficient $k$.

The smaller the value of $k$ the slower the decay in the spectrum.

In [6]:
gamma = 1e-5
delta = 1e-9

k = dl.Constant(1.0)
c = dl.Constant(0.1)

lmbda1, U1, Vm1, niter1 = solve(nx,ny, targets, gamma, delta,verbose=False)
k = dl.Constant(0.1)
lmbda2, U2, Vm2, niter2 = solve(nx,ny, targets, gamma, delta,verbose=False)
k = dl.Constant(0.01)
lmbda3, U3, Vm3, niter3 = solve(nx,ny, targets, gamma, delta,verbose=False)

print ("Number of Iterations: ", niter1, niter2, niter3)
plt.figure(figsize=(18,4))
nb.plot_eigenvalues(lmbda1, mytitle="Eigenvalues k=1.0", subplot_loc=131)
plt.ylim([1e-2,1e14])
nb.plot_eigenvalues(lmbda2, mytitle="Eigenvalues k=0.1", subplot_loc=132)
plt.ylim([1e-2,1e14])
nb.plot_eigenvalues(lmbda3, mytitle="Eigenvalues k=0.01", subplot_loc=133)
plt.ylim([1e-2,1e14])

nb.plot_eigenvectors(Vm1, U1, mytitle="k=1. Eigenvector", which=[0,1,5])
nb.plot_eigenvectors(Vm2, U2, mytitle="k=0.1 Eigenvector", which=[0,1,5])
nb.plot_eigenvectors(Vm3, U3, mytitle="k=0.01 Eigenvector", which=[0,1,5])

plt.show()

Number of Iterations:  38 129 775

7. Dependence on the number of observations

In [7]:
ntargets = [16, 64, 256]

gamma = 1e-5
delta = 1e-9

k = dl.Constant(0.1)
c = dl.Constant(0.1)

lmbda1, U1, Vm1, niter1 = solve(nx,ny, targets[0:ntargets[0],:], gamma, delta,verbose=False)
lmbda2, U2, Vm2, niter2 = solve(nx,ny, targets[0:ntargets[1],:], gamma, delta,verbose=False)
lmbda3, U3, Vm3, niter3 = solve(nx,ny, targets[0:ntargets[2],:], gamma, delta,verbose=False)

print ("Number of Iterations: ", niter1, niter2, niter3)
plt.figure(figsize=(18,4))
nb.plot_eigenvalues(lmbda1, mytitle="Eigenvalues ntargets={0}".format(ntargets[0]), subplot_loc=131)
plt.ylim([1e-6,1e12])
nb.plot_eigenvalues(lmbda2, mytitle="Eigenvalues ntargets={0}".format(ntargets[1]), subplot_loc=132)
plt.ylim([1e-6,1e12])
nb.plot_eigenvalues(lmbda3, mytitle="Eigenvalues ntargets={0}".format(ntargets[2]), subplot_loc=133)
plt.ylim([1e-6,1e12])

nb.plot_eigenvectors(Vm1, U1, mytitle="ntargets={0} Eigenvector".format(ntargets[0]), which=[0,1,5])
nb.plot_eigenvectors(Vm2, U2, mytitle="ntargets={0} Eigenvector".format(ntargets[1]), which=[0,1,5])
nb.plot_eigenvectors(Vm3, U3, mytitle="ntargets={0} Eigenvector".format(ntargets[2]), which=[0,1,5])

plt.show()

Number of Iterations:  31 94 129

Copyright © 2022, The University of Texas at Austin.

All Rights reserved. See file COPYRIGHT for details.

This file is part of cvips_labs, the teaching material for Computational and Variational Methods for Inverse Problems at The University of Texas at Austin. Please see https://hippylib.github.io/cvips_labs for more information and source code availability.

We would like to acknowledge the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program for providing cloud computing resources (Jetstream) for this course through allocation MTH230002. ACCESS is an advanced computing and data resource supported by the National Science Foundation and made possible through these lead institutions and their partners – Carnegie Mellon University; University of Colorado Boulder; University of Illinois at Urbana-Champaign; and State University of New York at Buffalo.