# An open-source framework for ExpFinder integrating $N$-gram Vector Space
Model and $\mu$CO-HITS

Hung Du<EMAIL_ADDRESS>Yong-Bin Kang<EMAIL_ADDRESS>Department of
Computer Science and Software Engineering, Swinburne University of Technology,
Australia  Department of Media and Communication, Swinburne University of
Technology, Australia

###### Abstract

Finding experts drives successful collaborations and high-quality product
development in academic and research domains. To contribute to the expert
finding research community, we have developed ExpFinder which is a novel
ensemble model for expert finding by integrating an $N$-gram vector space
model ($n$VSM) and a graph-based model ($\mu$CO-HITS). This paper provides
descriptions of ExpFinder’s architecture, key components, functionalities, and
illustrative examples. ExpFinder is an effective and competitive model for
expert finding, significantly outperforming a number of expert finding models
as presented in [1].

###### keywords:

ExpFinder , Expert finding , N-gram Vector Space Model , µCO-HITS , Expert
collaboration graph

††journal: Software Impacts

##  1 Introduction

Figure 1: The overview of ExpFinder

Identifying experts given a query topic, known as expert finding, is a crucial
task that accelerates rapid team formation for research innovations or
business growth. Existing expert finding models can be classified into three
categories such as vector space models (VSM) [2, 3], document language models
(DLM) [4, 5, 6], or graph-based models (GM) [7, 8, 9]. ExpFinder [1] is an
ensemble model for expert finding which integrates a novel $N$-gram VSM
($n$VSM) with a GM ($\mu$CO-HITS)-a variant of the generalised CO-HITS
algorithm [7].

As seen in Figure 1, ExpFinder has $n$VSM, a vector space model, as a key
component that estimates the weight of an expert and a document given a topic
by leveraging the Inverse Document Frequency (IDF) weighting [10] for $N$-gram
words (simply $N$-grams). Another key component in ExpFinder is $\mu$CO-HITS
that is used to reinforce the weights of experts and documents given a topic
in $n$VSM using an Expert Collaboration Graph (ECG) that is a certain form of
an expert social network. The output of ExpFinder is the reinforced weights of
experts given topics.

ExpFinder is designed and developed to improve the performance for expert
finding. In this paper, we highlight two main contributions to the expert
finding community. First, we provide a comprehensive implementation detail of
all steps taken in ExpFinder. It could also be used as an implementation
guideline for developing various DLM-, VSM- and GM-based expert finding
approaches. Second, we illustrate how ExpFinder works with a simple example,
thus researchers and practitioners can easily understand ExpFinder’s design
and implementation.

This paper is organised as follows. Section 2 describes ExpFinder’s
architecture and functionalities. Section 3 demonstrates the procedural steps
in ExpFinder. Section 4 provides the impact and conclusion of ExpFinder.

##  2 Functionality

ExpFinder is implemented in Python (version $\geq$ 3.6) with open-source
libraries such as pandas, NumPy, scikit-learn, SciPy, nltk, and networkx. In
this section, we present its architecture, key components, and their
functionalities. The architecture of ExpFinder is presented in Figure 2 that
consists of four key steps with the corresponding functions and their
functional dependencies:

Figure 2: The architecture and functional workflow of ExpFinder: blue labels
indicate module names of ExpFinder, and ‘Output Relation’ maps the functional
component to the corresponding processing step.

  1. 1.

Step 1 - Extract tokens and topics: Given an expertise source $\mathcal{D}$
(e.g., scientific publications) of experts $\mathcal{X}$, we extract expertise
topics by using tokenise_doc() in extractor.py. We assume that expertise
topics are represented in the forms of noun phrases. For each document
$d\in\mathcal{D}$, the function splits it into sentences. Then, for each
sentence, the function removes stopwords, assigns a part of speech (POS) to
each word, merges the inflected forms of a word (i.e., the lemmatisation
process, for example, ‘patients’ is lemmatised to ‘patient’), and extracts
single-word terms (called tokens) and topics with a given linguistic pattern.
In addition, we use a regular expression (regex) in Python to construct a
linguistic pattern based on POS that is further used for extracting four
different types of topics as shown in Figure 3.

Figure 3: The Python regular expression of a linguistic pattern for extracting
topics in a single document

Note that we use nltk for performing this process. The output of this step is
the list of the tokens and the list of topics for each document
$d\in\mathcal{D}$. The set of the all tokens is denoted as $\mathcal{W}$, and
the set of the all topics is denoted as $\mathcal{T}$.

  2. 2.

Step 2 - Estimate the weights of experts and documents given topics in $n$VSM:
The process includes four main steps with the corresponding functions in
generator.py:

    1. 2.1.

We use generate_tf() to estimate the term frequencies (TFs) of $\mathcal{W}$
in each document $d\in\mathcal{D}$. For this estimation, we use
CountVectoerizer in scikit-learn. The output of this function is the
$|\mathcal{D}|\times|\mathcal{W}|$ Document-Token matrix (DTM) where each
entry contains the TF of $w\in\mathcal{W}$ in $d$.

    2. 2.2.

We use generate_dp_matrix() to estimate the weights of documents given
$\mathcal{T}$ in $n$VSM [1]. The function estimates $n$TFIDF of each topic
$t\in\mathcal{T}$ by integrating the $n$TF weighting and the $n$IDF weighting.
Intuitively, $n$TF estimates the frequency of $t$ by averaging TFs of tokens
in $t$ where TF of each token is stored in DTM. In addition, $n$IDF [10] is
the $N$-gram IDF weighting method that estimates the log-IDF,
$\log\frac{|\mathcal{D}|\times df(t)+1}{df(w_{1}\land w_{2}\land\ldots\land
w_{n})^{2}+1}+1$, of $t$ where $w_{1},\ldots,w_{n}$ are $n$-constituent terms
in $t$. The output of this step is the $|\mathcal{D}|\times|\mathcal{T}|$
Document-Phrase matrix (DPM) where each entry contains the $n$TFIDF weight of
$t$ in $d$.

    3. 2.3.

Given $\mathcal{D}$, we use genenerate_ed_matrix() to generate the
$|\mathcal{X}|\times|\mathcal{D}|$ Expert-Document matrix (EDM) where each
entry shows a binary relationship between $x\in\mathcal{X}$ and $d$ (e.g., 1
indicates that $x$ has the authorship on $d$, and 0 otherwise).

    4. 2.4.

We use generate_pr_matrix() to estimate the weights of experts $\mathcal{X}$
and documents $\mathcal{D}$ given each topic $t\in\mathcal{T}$ in $n$VSM [1].
The weights of $\mathcal{X}$ are estimated by calculating matrix
multiplication of $\textbf{EDM}^{|\mathcal{X}|\times|\mathcal{D}|}$ and
$\textbf{DPM}^{|\mathcal{D}|\times|\mathcal{T}|}$ (e.g., ETopM =
numpy.matmul(EDM, DPM) in Python). The output is the
$|\mathcal{X}|\times|\mathcal{T}|$ Expert-Topic matrix (ETopM) where each
entry contains the topic-sensitive weight of $x$ given $t$. The weights of
$\mathcal{D}$ are represented by DPM. Now, we denote DPM as the
$|\mathcal{D}|\times|\mathcal{T}|$ Document-Topic matrix (DTopM) where each
entry shows the topic-sensitive weight of $d$ given $t$. It is worth noting
that DPM can be integrated with another factor (e.g., the average document
frequencies of $\mathcal{T}$) to obtain different weights for DTopM. However,
in our approach, we set $\textbf{DTopM}=\textbf{DPM}$.

  3. 3.

Step 3 - Construct ECG. We use generate_ecg() in generator.py to handle this
step. The function receives $\mathcal{D}$ and builds an ECG using DiGraph in
networkx to present a directed, weighted bipartite graph that has expert nodes
$V_{x}$ and document nodes $V_{d}$. The set of nodes in the graph is denoted
as $V$ such that $V=V_{x}\cup V_{d}$. A directed edge points from a document
node $v_{d}\in V_{d}$ to an expert node $v_{x}\in V_{x}$ if $x$ has published
$d$. In this step, we also use generate_ed_vector() in generator.py to
generate a $|V|\times 1$ Expert-Count vector ($\boldsymbol{c_{x}}$) and a
$|V|\times 1$ Document-Count ($\boldsymbol{c_{d}}$) vector based on ECG. These
vectors are used for the estimation of $\mu$CO-HITS in Step 4.

  4. 4.

Step 4 - Reinforce expert weights using $\mu$CO-HITS. We use run_expfinder()
in trainer.py to handle this step. The function receives ETopM, DTopM, ECG,
$\boldsymbol{c_{x}}$ and $\boldsymbol{c_{d}}$, generated in (Steps 2 and 3) as
parameters, and reinforces the estimation of expert weights given topics by
integrating $n$VSM and $\mu$CO-HITS [1]. For each $t\in\mathcal{T}$, we
perform the three steps:

    1. 4.1.

Generate the adjacency matrix of nodes and its transpose \- Given the ECG, we
use to_matrix() in networkx to generate the $|V|\times|V|$ adjacency matrix of
the graph M, and also construct its transpose matrix $\textbf{M}^{\top}$.
These matrices are required in the initialisation for running the $\mu$CO-HITS
algorithm.

    2. 4.2.

Normalise the weights of experts and documents given a topic \- We get topic-
sensitive weights of $\mathcal{X}$ and $\mathcal{D}$ given $t$ from ETopM and
DTopM, respectively. The output of this includes the $|\mathcal{X}|\times 1$
Expert-Topic ($\boldsymbol{\alpha_{x}}$) and $|\mathcal{D}|\times 1$ Document-
Topic ($\boldsymbol{\alpha_{d}}$) vectors where each entry shows the topic-
sensitive weight of an expert and a document given $t$, respectively. Then, we
normalise these vectors using L2 normalisation to scale their squares sum to 1
as the initialisation for running the $\mu$CO-HITS algorithm [11].

    3. 4.3.

Reinforce expert weights given a topic \- We integrate $n$VSM and $\mu$CO-HITS
through $k$ iterations to reinforce expert weights given $t$. $\mu$CO-HITS is
the extension of the CO-HITS algorithm [7] which contains two main properties
such as average authorities $\boldsymbol{a}$ and average hubs $\boldsymbol{h}$
which show importance of $\mathcal{X}$ and $\mathcal{D}$, respectively, based
on the ECG. In addition, these properties can be defined as [1]:

$\displaystyle\boldsymbol{a}(\mathcal{X};t)^{k}$
$\displaystyle=(1-\lambda_{x})\boldsymbol{a}(\mathcal{X};t)^{k-1}+\lambda_{x}\left(\frac{\textbf{M}^{\top}\boldsymbol{\cdot}\boldsymbol{h}(\mathcal{D};t)^{k-1}}{\boldsymbol{c_{d}}}\right)$
(1) $\displaystyle\boldsymbol{h}(\mathcal{D};t)^{k}$
$\displaystyle=(1-\lambda_{d})\boldsymbol{h}(\mathcal{D};t)^{k-1}+\lambda_{d}\left(\frac{\textbf{M}\boldsymbol{\cdot}\boldsymbol{a}(\mathcal{X};t)^{k}}{\boldsymbol{c_{x}}}\right)$
(2)

where

       * •

$\boldsymbol{a}(\mathcal{X};t)^{k}$ and $\boldsymbol{h}(\mathcal{D};t)^{k}$
are $|V|\times 1$ vectors which contain the reinforced expert weights and
reinforced document weights, respectively, given $t$ at $k^{\text{th}}$
iteration. As the initial weights of these vectors, we use the topic-sensitive
weights of experts and documents estimated in $n$VSM. Thus,
$\boldsymbol{a}(\mathcal{X};t)^{0}=\boldsymbol{\alpha_{x}}$ and
$\boldsymbol{h}(\mathcal{D};t)^{0}=\boldsymbol{\alpha_{d}}$. By doing so, we
integrate $n$VSM with $\mu$CO-HITS. Note that
$\boldsymbol{a}(\mathcal{X};t)^{0}$ is a $|\mathcal{X}|\times 1$ vector, and
$\boldsymbol{h}(\mathcal{D};t)^{0}$ is a $|\mathcal{D}|\times 1$ vector.
However, for easily implementing the HITS algorithm, we have transformed the
dimension of these vectors into $|V|\times 1$ vectors where additional entries
hold the value of 0.

       * •

$\lambda_{x}\in[0,1]$ and $\lambda_{d}\in[0,1]$ are parameters for expert and
document, respectively. These are used to control the impact of topic-
sensitive weights on $\boldsymbol{a}$ and $\boldsymbol{h}$, respectively.
Assigning lower values indicates the higher impact of topic-sensitive weights
on $\boldsymbol{a}$ and $\boldsymbol{h}$.

       * •

$\left(\frac{\textbf{M}^{\top}\boldsymbol{\cdot}\boldsymbol{h}(\mathcal{D};t)^{k-1}}{\boldsymbol{c_{d}}}\right)$
is the calculation for the average authorities. The numerator performs matrix
multiplication between the $|V|\times|V|$ adjacency matrix $\textbf{M}^{\top}$
and the $|V|\times 1$ $\boldsymbol{h}$. The denominator is a $|V|\times 1$
counted vector $\boldsymbol{c_{d}}$ generated in Step 3. To calculate this in
Python, we simply apply numpy.matmul($\textbf{M}^{\top}$,
$\boldsymbol{h}(\mathcal{D};t)^{k-1}$)/$\boldsymbol{c_{d}}$.

       * •

$\left(\frac{\textbf{M}\boldsymbol{\cdot}\boldsymbol{a}(\mathcal{X};t)^{k}}{\boldsymbol{c_{x}}}\right)$
is the calculation for the average hubs. The numerator performs matrix
multiplication between the $|V|\times|V|$ adjacency matrix M and the
$|V|\times 1$ $\boldsymbol{a}$. The denominator is a $|V|\times 1$ counted
vector $\boldsymbol{c_{x}}$ generated in Step 3. To calculate this in Python,
we simply apply numpy.matmul(M,
$\boldsymbol{a}(\mathcal{X};t)^{k}$)/$\boldsymbol{c_{x}}$.

After computing $\boldsymbol{a}$ and $\boldsymbol{h}$ at $k^{\text{th}}$
iteration, we apply L2 normalisation to both $\boldsymbol{a}$ and
$\boldsymbol{h}$. We use the obtained $\boldsymbol{a}(\mathcal{X};t)$ after
the final iteration to construct the $|\mathcal{X}|\times|\mathcal{T}|$
Expert-Topic matrix (RETopM) where each entry contains the reinforced weight
of $x$ given $t$.

##  3 Illustrative examples

In this section, we illustrate how ExpFinder works. The input data111The
example data are also provided in our Github repository. includes three
experts (i.e., $x_{1}$, $x_{2}$ and $x_{3}$) and three documents (i.e.,
$d_{1}$, $d_{2}$ and $d_{3}$) as shown in Table 1. Figure 4 presents the
output examples of the steps in ExpFinder:

Docs |  Experts |  Text
---|---|---
$d_{1}$ |  $x_{1}$, $x_{2}$ |  A prerequisite for using electronic health records (EHR) data within learning health-care system is an infrastructure that enables access to EHR data longitudinally for health-care analytics and real time for knowledge delivery. Herein, we share our institutional implementation of a big data-empowered clinical natural language processing (NLP) infrastructure, which not only enables healthcare analytics but also has real-time NLP processing capability.
$d_{2}$ |  $x_{1}$, $x_{3}$ |  Word embedding, where semantic and syntactic features are captured from unlabeled text data, is a basic procedure in Natural Language Processing (NLP). In this paper, we first introduce the motivation and background of word embedding and its related language models.
$d_{3}$ |  $x_{2}$ |  Structural health monitoring at local and global levels using computer vision technologies has gained much attention in the structural health monitoring community in research and practice. Due to the computer vision technology application advantages such as non-contact, long distance, rapid, low cost and labor, and low interference to the daily operation of structures, it is promising to consider computer vision structural health monitoring as a complement to the conventional structural health monitoring. This article presents a general overview of the concepts, approaches, and real-life practice of computer vision structural health monitoring along with some relevant literature that is rapidly accumulating.

Table 1: The document dataset $\mathcal{D}$ used in the example: extracted
phrases are highlighted in yellow, and extracted tokens are in bold. Figure 4:
Illustrative examples for ExpFinder: the blue labels indicate module names of
ExpFinder.

  1. 1.

Step 1 - Extract tokens and topics: Given $\mathcal{D}$, we extract tokens
$\mathcal{W}$ and topics $\mathcal{T}$. In this step, we set a maximum length
of phrase to be 3 such that we only obtain phrases that have less than or
equal to 3 tokens. Additionally, we use the linguistic pattern presented in
Section 2. The output of this step contains the set of 50 unique topics
$\mathcal{T}$ and the set of 85 unique tokens $\mathcal{W}$. For example,
extracted topics in $d_{1}$ include some single-token topics (e.g.,
prerequisite and capability) and some multi-token topics (e.g., real-time nlp
processing and electronic health record).

  2. 2.

Step 2 - Estimate the weights of experts and documents given topics \- Given
$\mathcal{T}$ and $\mathcal{W}$, we generate three main matrices (i.e., EDM,
DTopM and ETopM) that will also be used in Step 4. To do this, we perform the
following:

     * •

Given $\mathcal{W}$, we generate $\textbf{DTM}^{3\times 85}$ where each entry
shows the TF of a token $w\in\mathcal{W}$ in a document $d\in\mathcal{D}$. For
example, the $3\times 1$ vector of healthcare,
$\textbf{DTM}_{*,\texttt{healthcare}}$, is $(1,0,0)$ which shows it occurs
only in $d_{1}$ (see also $\mathcal{D}$ in Table 1). As another example, we
obtain $\textbf{DTM}_{*,\texttt{analytics}}=(2,0,0)$ which denotes that
analytics appears twice in $d_{1}$.

     * •

Given $\mathcal{T}$ and $\textbf{DTM}^{3\times 85}$, we generate
$\textbf{DPM}^{3\times 50}$ where each entry contains the weight of a phrase
$t\in\mathcal{T}$ for a document calculated in $n$VSM. For example, suppose
that health analytics is denoted as $t_{1}$, we then calculate $n$TF of
$t_{1}$ in $d_{1}$ as:

$\displaystyle
n\text{TF}(t_{1},d_{1})=\frac{\textbf{DTM}_{1,\texttt{healthcare}}+\textbf{DTM}_{1,\texttt{analytics}}}{|t_{1}|}\
=\frac{(1+2)}{2}=1.5$

where $|t_{1}|$ is a number of tokens in $t_{1}$. Then, we calculate the
$N$-gram IDF of $t_{1}$ as:

$\displaystyle n\text{IDF}(t_{1})$ $\displaystyle=\log{\frac{|D|\cdot
df(t_{1})+1}{df(\textbf{DTM}_{*,\texttt{healthcare}}\land\textbf{DTM}_{*,\texttt{analytics}})^{2}+1}}+1$
$\displaystyle=\log{\frac{3\times 1+1}{df((1,0,0)\land(2,0,0))^{2}+1}}+1$
$\displaystyle=\log{\frac{4}{1^{2}+1}}+1=1.693.$

Here, $\land$ is implemented in NumPy. Finally, we multiply
$n\text{TF}(t_{1},d_{1})$ with $n\text{IDF}(t_{1})$ to obtain $n$TFIDF of
$d_{1}$ given $t_{1}$ as:
$\textbf{DPM}_{1,t_{1}}=n\text{TFIDF}(t_{1},d_{1})=n\text{TF}(t_{1},d_{1})\times
n\text{IDF}(t_{1})=1.5\times 1.693=2.540$.

     * •

Given $\mathcal{D}$, we generate $\textbf{EDM}^{3\times 3}$ where each entry
shows the authorship of an expert on a document. For example, $x_{1}$ is an
author of $d_{1}$, and hence, the entry between $x_{1}$ and $d_{1}$
($\textbf{EDM}_{1,1}$) equals 1. Also, $\textbf{EDM}_{1,3}=0$ shows that
$x_{1}$ is not an author of $d_{3}$ (See Table 1).

     * •

Given $\textbf{EDM}^{3\times 3}$ and $\textbf{DPM}^{3\times 50}$, we generate
$\textbf{ETopM}^{3\times 50}$ where each entry contains $n$TFIDF weight of an
expert given a topic. As we explained in Section 2, we assume that
$\textbf{DTopM}=\textbf{DPM}$. Now, we demonstrate the calculation for the
weights of experts $\mathcal{X}$ given $t_{1}$ ($\textbf{ETopM}_{*,t_{1}}$) in
$n$VSM as:

$\displaystyle\textbf{ETopM}_{*,t_{1}}$ $\displaystyle=\textbf{EDM}^{3\times
3}\boldsymbol{\cdot}\textbf{DTopM}_{*,t_{1}}$
$\displaystyle=\begin{bmatrix}1&1&0\\\ 1&0&1\\\
0&1&0\end{bmatrix}\boldsymbol{\cdot}(2.540,0,0)=(2.540,2.540,0)$

Note that $\textbf{DTopM}^{3\times 3}$ and $\textbf{ETopM}^{3\times 3}$ are
only used for the visualisation purpose. We use $\textbf{DTopM}^{3\times 50}$
and $\textbf{ETopM}^{3\times 50}$ for the estimation in Step 4.

  3. 3.

Step 3 - Construct ECG: Given $\mathcal{D}$, we generate an ECG which has
three expert nodes and three document nodes, as shown in Figure 4. The graph
is also used to generate $3\times 1$ vectors (i.e., $\boldsymbol{c_{x}}$ and
$\boldsymbol{c_{d}}$) that are used for the estimation of $\mu$CO-HITS in Step
4. For example, $\boldsymbol{c_{d}}_{1}=2$ indicates there are two documents
(i.e., $d_{1}$ and $d_{2}$) pointing to $x_{1}$. Similarly,
$\boldsymbol{c_{x}}_{3}=1$ indicates that there is one expert (i.e., $x_{2}$)
who has authorship on $d_{3}$.

  4. 4.

Step 4 - Reinforce expert weights using $\mu$CO-HITS: We use run_expfinder()
in trainer.py to reinforce expert weights given topics $\mathcal{T}$. The
function receives $\textbf{DTopM}^{3\times 50}$, $\textbf{ETopM}^{3\times
50}$, ECG, $\boldsymbol{c_{x}}$ and $\boldsymbol{c_{d}}$, generated in (Steps
2 and 3) as parameters, and generate the $3\times 50$ Expert-Topic matrix
where each entry shows the reinforced weight of an expert given a topic. Now,
we illustrate the estimation for the reinforced weight of $\mathcal{X}$ given
$t_{1}$ as:

     * •

Given 6 nodes in an ECG, we generate the adjacency matrix $\textbf{M}^{6\times
6}$ and its transpose matrix $\textbf{M}^{\top}$ as:

$\displaystyle\textbf{M}=\begin{bmatrix}0&1&1&0&0&0\\\ 0&0&0&0&0&0\\\
0&0&0&0&0&0\\\ 0&1&0&0&1&0\\\ 0&0&0&0&0&0\\\
0&0&1&0&0&0\end{bmatrix},\textbf{M}^{\top}=\begin{bmatrix}0&0&0&0&0&0\\\
1&0&0&1&0&0\\\ 1&0&0&0&0&1\\\ 0&0&0&0&0&0\\\ 0&0&0&1&0&0\\\
0&0&0&0&0&0\end{bmatrix}$

where rows and columns are labeled with the sequence $\boldsymbol{s}$ (i.e.,
$\boldsymbol{s}=(d_{1},x_{1},x_{2},d_{2},x_{3},d_{3})$).

     * •

We apply L2 normalisation for the $6\times 1$ Expert-Topic
($\boldsymbol{\alpha_{x}}$) and the $6\times 1$ Document-Topic
($\boldsymbol{\alpha_{d}}$) vectors. The output of each vector is as:

$\displaystyle\boldsymbol{\alpha_{x}}$
$\displaystyle=\texttt{L2-normalize}(\textbf{ETopM}_{*,t_{1}})=(0,0.707,0.707,0,0,0)$
$\displaystyle\boldsymbol{\alpha_{d}}$
$\displaystyle=\texttt{L2-normalize}(\textbf{DTopM}_{*,t_{1}})=(1,0,0,0,0,0)$

     * •

We reinforce expert weights given $t_{1}$ in 5 iterations with $\lambda_{x}=1$
and $\lambda_{d}=0.7$. Here, we demonstrate the calculation of average
authorities $\boldsymbol{a}$ and average hubs $\boldsymbol{h}$ at the first
iteration ($k=1$):

$\displaystyle\boldsymbol{a}(\mathcal{X};t_{1})^{1}$
$\displaystyle=(1-\lambda_{x})\boldsymbol{a}(\mathcal{X};t_{1})^{0}+\lambda_{x}\left(\frac{\textbf{M}^{\top}\boldsymbol{\cdot}\boldsymbol{h}(\mathcal{D};t_{1})^{0}}{\boldsymbol{c_{d}}}\right)$
$\displaystyle=0\cdot(0,0.707,0.707,0,0,0)+1.0\cdot\left(\frac{(0,2,2,0,1,0)}{(2,1,1,2,1,1)}\right)$
$\displaystyle=(0,2,2,0,1,0)$
$\displaystyle\boldsymbol{h}(\mathcal{D};t_{1})^{1}$
$\displaystyle=(1-\lambda_{d})\boldsymbol{h}(\mathcal{D};t_{1})^{0}+\lambda_{d}\left(\frac{\textbf{M}\boldsymbol{\cdot}\boldsymbol{a}(\mathcal{X};t_{1})^{1}}{\boldsymbol{c_{x}}}\right)$
$\displaystyle=0.3\cdot(1,0,0,0,0,0)+0.7\cdot\left(\frac{(4,0,0,3,0,2)}{(1,2,2,1,1,1)}\right)$
$\displaystyle=(3.1,0,0,2.1,0,1.4)$

where $\boldsymbol{a}(\mathcal{X};t_{1})^{1}$ and
$\boldsymbol{h}(\mathcal{D};t_{1})^{1}$ are $6\times 1$ vectors. At the end of
the iteration, we normalise these vectors by applying the L2 normalisation
technique as:

$\displaystyle\boldsymbol{a}(\mathcal{X};t_{1})^{1}$
$\displaystyle=\texttt{L2-normalize}(\boldsymbol{a}(\mathcal{X};t_{1})^{1})=(0,0.667,0.667,0,0.333,0)$
$\displaystyle\boldsymbol{h}(\mathcal{D};t_{1})^{1}$
$\displaystyle=\texttt{L2-normalize}(\boldsymbol{h}(\mathcal{D};t_{1})^{1})=(0.776,0,0,0.525,0,0.35)$

After 5 iterations, we obtain
$\boldsymbol{a}(\mathcal{X};t_{1})^{5}=(0,0.577,0.595,0,0.56,0)$ whose labels
are presented by $\boldsymbol{s}$, and hence, we use $x_{1}$, $x_{2}$ and
$x_{3}$ as indexes for obtaining a $3\times 1$ vector (i.e.,
$\textbf{RETopM}_{*,t_{1}}=(0.577,0.595,0.56)$).

The output is $\textbf{RETopM}^{3\times 50}$. If we use $t_{1}$ and the other
two topics (i.e., natural language processing and vision technology, denoted
as $t_{2}$ and $t_{3}$, respectively), we can generate
$\textbf{RETopM}^{3\times 3}$ in Figure 4. This matrix can be used for two
major tasks (1) finding the most expertise query for each expert (also known
as expert profiling); and (2) finding the best expert for a given query (also
known as expert finding).

##  4 Impact and Conclusion

With the growth of expertise digital sources, expert finding is a crucial task
that has significantly helped people to seek the services and guidance of an
expert [12]. ExpFinder is an ensemble model for expert finding that integrates
$n$VSM with $\mu$CO-HITS to enhance the capability for expert finding over
existing DLM, VSM and GM approaches. To our best knowledge, ExpFinder is the
first attempt to provide the implementation of $n$VSM and $\mu$CO-HITS for
expert finding.

The implementation of ExpFinder also provides functionalities that can be
potentially useful for implementing other expert finding models. For example,
our tokenisation module for extracting noun phrases using a linguistic pattern
based on a part of speech (POS) can be easily customised based on researchers’
purposes. The modules for building the presented Expert-Document matrix (EDM),
Expert-Topic matrix (ETopM), and Document-Topic matrix (DTopM) can be usefully
leveraged to represent relationships between experts and documents, experts
and topics, and documents and topics. These relationships can be used to
represent a collective information among experts, documents and topics and
used to implement other graph-based expert finding models such as an author-
document-topic (ADT) graph [8] and an expert-expert graph via topics.

We highlight that ExpFinder is a state-of-the-art model, substantially
outperforming the following widely known and latest models for expert finding:
document language models [4], probabilistic-based expert finding model [6]),
graph-based models [7, 8, 9]. Thus, the ones who want to extend ExpFinder can
harness our implementation for further improvement of ExpFinder.

We presented the architecture and implementation detail of ExpFinder with an
illustrative example. This would help researchers and practitioners to better
understand how ExpFinder is designed and implemented with its core
functionalities.

##  5 Declaration of Competing Interest

The authors declare that they have no known competing financial interests or
personal relationships that could have appeared to influence the work reported
in this paper.

## References

  * [1] Y.-B. Kang, H. Du, A. R. M. Forkan, P. P. Jayaraman, A. Aryani, T. Sellis, Expfinder: An ensemble expert finding model integrating $n$-gram vector space model and $\mu$co-hits (2021).  arXiv:2101.06821.
  * [2] F. Riahi, Z. Zolaktaf, M. Shafiei, E. Milios, Finding expert users in community question answering, in: Proceedings of the 21st International Conference on World Wide Web, 2012, pp. 791–798.
  * [3] C. T. Chuang, K. H. Yang, Y. L. Lin, J. H. Wang, Combining query terms extension and weight correlative for expert finding, in: 2014 IEEE/WIC/ACM International Joint Conferences on Web Intelligence (WI) and Intelligent Agent Technologies (IAT), Vol. 1, IEEE, 2014, pp. 323–326.
  * [4] K. Balog, L. Azzopardi, M. de Rijke, A language modeling framework for expert finding, Information Processing & Management 45 (1) (2009) 1 – 19.
  * [5] B. Wang, X. Chen, H. Mamitsuka, S. Zhu, Bmexpert Mining medline for finding experts in biomedical domains based on language model, IEEE/ACM Trans. Comput. Biol. Bioinformatics 12 (6) (2015) 1286–1294.
  * [6] P. Cifariello, P. Ferragina, M. Ponza, WISER: A semantic approach for expert finding in academia based on entity linking, Information Systems 82 (2019) 1 – 16.
  * [7] H. Deng, M. R. Lyu, I. King, A Generalized CO-HITS Algorithm and Its Application to Bipartite Graphs, in: Proceedings of the 15th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, 2009, p. 239–248.
  * [8] S. D. Gollapalli, P. Mitra, C. L. Giles, Ranking experts using author-document-topic graphs, in: Proceedings of the 13th ACM/IEEE-CS Joint Conference on Digital Libraries, JCDL ’13, 2013, pp. 87–96.
  * [9] D. Schall, A Social Network-Based Recommender Systems, Springer, 2015.
  * [10] M. Shirakawa, T. Hara, S. Nishio, IDF for Word N-Grams, ACM Trans. Inf. Syst. 36 (1).
  * [11] J. M. Kleinberg, Authoritative sources in a hyperlinked environment, Journal of the ACM (JACM) 46 (5) (1999) 604–632.
  * [12] R. Gonçalves, C. F. Dorneles, Automated expertise retrieval: A taxonomy-based survey and open issues, ACM Computing Surveys (CSUR) 52 (5) (2019) 1–30.

##  6 Current code version

Ancillary data table required for subversion of the codebase. Kindly replace
examples in right column with the correct information about your current code,
and leave the left column as it is.

Nr. |  Code metadata description |
---|---|---
C1 |  Current code version |  v1.0
C2 |  Permanent link to code/repository used for this code version |  https://github.com/Yongbinkang/ExpFinder
C3 |  Permanent link to Reproducible Capsule |  https://doi.org/10.24433/CO.0133456.v1
C4 |  Legal Code License |  MIT License (MIT)
C5 |  Code versioning system used |  git
C6 |  Software code languages, tools, and services used |  Python
C7 |  Compilation requirements, operating environments & dependencies |  Python environment version 3.6 or above, pandas, networkx, NumPy, scikit-learn, nltk, SciPy, Torch, Transformers, SciBERT
C8 |  Link to developer documentation/manual |  https://github.com/Yongbinkang/ExpFinder/blob/main/README.md
C9 |  Support email for questions |<EMAIL_ADDRESS><EMAIL_ADDRESS>
Table 2: Code metadata


# Inverse Problems with Diffusion Models: A MAP Estimation Perspective

Sai bharath chandra Gutha
KTH, Sweden
<EMAIL_ADDRESS>Hossein Azizpour
KTH, Sweden
<EMAIL_ADDRESS>Ricardo Vinuesa
KTH, Sweden
<EMAIL_ADDRESS>
###### Abstract

Inverse problems have many applications in science and engineering. In
Computer vision, several image restoration tasks such as inpainting,
deblurring, and super-resolution can be formally modeled as inverse problems.
Recently, methods have been developed for solving inverse problems that only
leverage a pre-trained unconditional diffusion model and do not require
additional task-specific training. In such methods, however, the inherent
intractability of determining the conditional score function during the
reverse diffusion process poses a real challenge, leaving the methods to
settle with an approximation instead, which affects their performance in
practice. Here, we propose a MAP estimation framework to model the reverse
conditional generation process of a continuous time diffusion model as an
optimization process of the underlying MAP objective, whose gradient term is
tractable. In theory, the proposed framework can be applied to solve general
inverse problems using gradient-based optimization methods. However, given the
highly non-convex nature of the loss objective, finding a perfect gradient-
based optimization algorithm can be quite challenging, nevertheless, our
framework offers several potential research directions. We use our proposed
formulation and develop empirically effective algorithms for solving noiseless
and noisy image inpainting tasks. We validate our proposed algorithms with
extensive experiments across diverse mask settings.

##  1 Introduction

Inverse problems are ubiquitous in science and engineering with a wide range
of downstream applications [27, 2]. In Computer vision, several image
restoration tasks such as inpainting, deblurring, super-resolution, and more,
can be formally modeled as inverse problems. In an inverse problem,
characterized by Eq. 1, $y\in\mathbb{R}^{m}$ is a (potentially noisy)
observation of the original data $x\in\mathbb{R}^{n}$, and $\eta$ is a random
variable denoting i.i.d noise, typically assumed to be Gaussian with a known
variance i.e $\eta\sim\mathcal{N}(0,\sigma^{2}_{y}\mathbb{I})$, and the task
is to infer the original data $x$ given the observation $y$. The function
$\mathcal{A}:\mathbb{R}^{n}\to\mathbb{R}^{m}$ is known as the forward
operator, and typically $n\gg m$, indicating that the observation
$y\in\mathbb{R}^{m}$ corresponds to a severely degraded signal, from which one
needs to recover the original signal $x\in\mathbb{R}^{n}$, which makes the
task highly challenging. For linear inverse problems, $\mathcal{A}$ denotes a
linear mapping and can be substituted with a matrix $H\in\mathbb{R}^{m\times
n}$.

$y=\mathcal{A}(x)+\eta$ (1)

Several conventional approaches for solving inverse problems exist [1]. These
include approaches based on functional-analytic, probabilistic (bayesian),
data-driven methods, and more. Recently, Deep Learning (DL) based methods have
been applied to solve inverse problems and have shown great success. In a
bayesian framework, solving an inverse problem naturally corresponds to
estimating the posterior $P(x|y)$. Typical DL-based approaches for solving
inverse problems fall into two categories. 1. Methods that directly learn the
posterior $P(x|y)$ via conditional generative models [13, 18], and 2. Methods
that learn $P(x)$ via an unconditional generative model and use it to infer
$P(x|y)$ [20, 5, 28, 11]. Methods of the former category require task-specific
training, i.e. training with a dataset of pairs $(x,y)$, where the degradation
$y$ is computed using $x$ and a task-specific forward operator $\mathcal{A}$.
This limits the out-of-the-box applicability of the model to a different task
(different forward operator). On the contrary, methods of the latter category
train an unconditional generative model to learn $P(x)$, and this training is
task-independent since it only needs a dataset of original data samples $x$.
These methods then use the trained model for $P(x)$ and since $P(y|x)$ is
tractable ( i.e from Eq. 1,
$P(y|x)=\mathcal{N}(\mathcal{A}(x),\sigma^{2}_{y}\mathbb{I})$ ), utilizing the
bayes rule, they infer the posterior $P(x|y)\propto P(y|x)P(x)$.

Several choices for Deep Generative Models (DGMs) exist, each with its
advantages and disadvantages. There have been approaches using Generative
Adversarial Networks (GAN) [8] and Normalizing Flow (NF) [16] based DGMs for
solving inverse problems, with more recent methods focusing on Diffusion
models [19, 9, 25], owing to their state-of-the-art performance in several
vision-based generative tasks. In this work, we focus on methods that use a
pre-trained unconditional diffusion model as the prior $P(x)$ and infer the
posterior $P(x|y)$ for solving inverse problems. Sec. 2 provides some
background on diffusion models and related works that use unconditional
diffusion models for solving inverse problems and their inherent limitations.
Later in Sec. 3, we propose our Maximum A Posteriori (MAP) estimation
framework for continuous-time diffusion models and discuss the practical
implementation. In Sec. 4, use our proposed framework to develop empirically
effective algorithms for image inpainting, and in Sec. 5, we validate our
proposed method with extensive experiments on image inpainting across diverse
mask settings. Finally, in Secs. 6 and 7, we discuss concurrent works and
conclude our findings.

##  2 Background

Given a dataset of samples $\mathcal{D}={\\{x_{i}\\}}_{i=1}^{N}$, where each
$x_{i}$ is an i.i.d sample drawn from an unknown data distribution
$P_{data}(x)$, a generative model learns to approximate $P_{data}$ from the
samples in $\mathcal{D}$.

###  2.1 Diffusion Models

Diffusion models [19, 9, 25] are a recent family of generative models. The
methodology involves simulating a stochastic process $\\{x(t)\\}_{t=0}^{T}$
described by an SDE (Stochastic Differential Equation) such as Eq. 2, where
$t\in[0,T]$ is a continuous time variable, $x(0)\sim P_{0}=P_{data}$ is the
data distribution for which we have a dataset $\mathcal{D}$ of samples,
$x(T)\sim P_{T}$ is a tractable prior distribution. The functions
$f(.,t):\mathbb{R}^{n}\to\mathbb{R}^{n}$ and $g(.):\mathbb{R}\to\mathbb{R}$
are called the drift and diffusion coefficients of $x(t)$ respectively, and
$dw$ denotes the standard Wiener process. Typically $f$ and $g$ are chosen in
a way that yields a tractable prior $P_{T}$ which contains no information
about $P_{0}$ (i.e. $P_{data}$).

$dx=f(x,t)dt+g(t)dw$ (2)

Eq. 2 also describes the ”forward process”, in which, starting from an
initially clean data sample, a small amount of noise is progressively added at
each step until it turns into a noisy sample of the prior distribution
$P_{T}$. The backward/reverse process which transforms a noisy sample of
$P_{T}$ into a clean sample of the data distribution is described by the
corresponding reverse-SDE in Eq. 3

$dx=[f(x,t)-g(t)^{2}\nabla_{x}\log P_{t}(x)]dt+g(t)d\bar{w}$ (3)

$d\bar{w}$ denotes the standard Wiener process, and $t$ goes backward from $T$
to $0$ with $dt$ being an infinitesimal negative time step. The term
$\nabla_{x}\log P_{t}(x)$ is called the score function of the marginal
distribution $P_{t}(x)$. If we know this score function for each marginal
distribution i.e. for all $t$, then one could, in theory, simulate the
reverse-SDE in Eq. 3 and generate samples from the data distribution. In
general, the score function is not analytically tractable and is hard to
estimate, however, one could train a time-indexed neural network model to
learn the score function via score matching techniques [26, 24]. The trained
score model $S_{\theta}(x,t)$ can then be substituted in place of
$\nabla_{x_{t}}\log P_{t}(x)$ in Eq. 3, and the reverse-SDE can be solved
using traditional SDE solvers [25].

$dx=[f(x,t)-\frac{1}{2}g(t)^{2}\nabla_{x}\log P_{t}(x)]dt$ (4)

For the stochastic process described by the SDE in Eq. 2, there exists a
corresponding ODE (Ordinary Differential Equation) shown in Eq. 4 describing a
deterministic process whose trajectories share the same marginal probability
densities $\\{P(x_{t})\\}_{t=0}^{T}$ as those simulated by the SDE. This is
referred to as the probability flow ODE in the literature [25]. So
equivalently, one could also use ODE solvers to solve Eq. 4 in reverse time
from $t=T$ until $0$ to generate samples from the data distribution.

From hereon, unless mentioned otherwise, we assume the default choice for
drift and diffusion coefficients as $f(x,t)=0$ and
$g(t)=\sqrt{\frac{d\sigma^{2}(t)}{dt}}$, where $\sigma(t)$ is a monotonically
increasing noise schedule from $t=0$ to $T$, with $\sigma(T)$ being very high.
This choice of $f$ and $g$ results in a closed form perturbation kernel
$P(x(t)|x(0))=\mathcal{N}(x(0),\sigma^{2}(t)-\sigma^{2}(0))$ and
$P(x(T))\approx\mathcal{N}(0,\sigma^{2}(T))$. There can be multiple design
choices for the noise schedule $\sigma(t)$, resulting in various formulations
of diffusion models [10].

###  2.2 Solving Inverse problems with Diffusion models

As described in Sec. 1, solving an inverse problem entails estimation of (or
sampling from) the posterior $P(x|y)$, where $y$ is the noisy degradation of
$x$ (Eq. 1). In the context of solving inverse problems using diffusion
models, sampling from the posterior $P(x|y)$ involves conditioning the reverse
diffusion process on $y$ which translates to solving the modified reverse-SDE
in Eq. 5.

$dx=[f(x,t)-g(t)^{2}\nabla_{x}\log P_{t}(x|y)]dt+g(t)d\bar{w}$ (5)

Similar to methods of the first category (Sec. 1), which directly learn the
posterior $P(x|y)$ as part of their training, it is possible to train a
conditional diffusion model that learns the conditional score function
directly. More specifically, one could learn $S_{\theta}(x,y,t)$ using
conditional score matching objectives in place of the usual unconditional
score function $S_{\theta}(x,t)$ in Sec. 2.1. In this paper, we focus on
methods of the second category, which only leverage an unconditional diffusion
model for $P(x)$, to infer $P(x|y)$. We modify the notations to denote
$P_{t}(x)$ with $P(x_{t})$ and $x(t)$ with $x_{t}$ for convenience.

$\nabla_{x_{t}}\log P(x_{t}|y)=\nabla_{x_{t}}\log P(x_{t})+\nabla_{x_{t}}\log
P(y|x_{t})$

(6)

$P(y|x_{t})=\int_{x_{0}}P(y|x_{0})P(x_{0}|x_{t})dx_{0}$

(7)

Solving Eq. 5 involves estimating the conditional score function
$\nabla_{x_{t}}\log P(x_{t}|y)$. The pre-trained unconditional diffusion model
can be used to estimate $\nabla_{x_{t}}\log P(x_{t})$, however, the term
$\nabla_{x_{t}}\log P(y|x_{t})$ becomes intractable, making it hard to
estimate the conditional score (Eqs. 6 and 7). At its core, the intractability
of $\nabla_{x_{t}}\log P(y|x_{t})$ arises from the fact that $P(x_{0}|x_{t})$
is intractable [20].

###  2.3 Related works

PGDM [20] approximates $P(x_{0}|x_{t})$ with a Gaussian distribution having
mean $\hat{x_{t}}$ and variance ${r_{t}}^{2}$, where
$\hat{x_{t}}=\mathbb{E}(x_{0}|x_{t})=x_{t}+\sigma^{2}(t)\nabla_{x_{t}}\log
P(x_{t})$ (using Tweedie’s formula). The standard deviation $r_{t}$ is a
hyperparameter, chosen proportionally to $\sigma(t)$. DPS [5] approximates
$P(y|x_{t})$ with the point estimate $P(y|x_{0}=\hat{x_{t}})$ and has an
almost similar formulation as PGDM, though the motivation is slightly
different. [3] also uses Gaussian approximation but further replace $r_{t}$
with the covariance matrix
$Cov[x_{0}|x_{t}]=\sigma^{2}(t)\frac{\partial\hat{x_{t}}}{\partial x_{t}}$.
Computing this matrix is expensive in practice, so they resort to diagonal and
row-sum approximations of the matrix instead. [15] proposes to find an optimal
covariance matrix using learned covariances from the diffusion model. All
those methods, however, still rely on making simplified approximations for
$P(x_{0}|x_{t})$. This limits their performance in practice, given the
complicated and multimodal nature of the true data distribution. Other lines
of work [11, 29, 7, 6] try to circumvent this term by projecting the
intermediate $x_{t}$ onto the measurement subspace using heuristic
approximations. Some other works [6, 31, 4, 28] pose the problem as a MAP
optimization with data and prior terms. Note that our proposed method is quite
similar to the concurrent works [4, 28] from a practical perspective. We
discuss the similarities and differences later in Sec. 6.

##  3 Our Methodology

###  3.1 Background: Consistency Models

Consider the probability flow ODE described in Eq. 4. The solution
trajectories of this ODE are smooth, and maps the samples on the data manifold
to pure noise. In [22], a consistency model is defined as the function that
map any point on the PF ODE trajectory to its corresponding origin (initial
point on the data manifold). There exist efficient methodologies [21] to train
these consistency models in practice. Please refer to [22] for a detailed
description of consistency models.

###  3.2 Proposed MAP estimation framework

$x^{*}_{0}=\arg\max_{x_{0}}\ \log P(x_{0}|y)$ (8)

Eq. 8 refers to the usual MAP formulation for solving an inverse problem.
Finding an optimal $x^{*}_{0}$ using gradient ascent involves the update step
in Eq. 9.

$x^{k+1}_{0}=x^{k}_{0}+\lambda*\nabla_{x_{0}}\log P(x_{0}|y)$ (9)

The update step requires computing the gradient term $\nabla_{x_{0}}\log
P(x_{0}|y)=\nabla_{x_{0}}\log P(y|x_{0})+\nabla_{x_{0}}\log P(x_{0})$. The
former term is tractable since $P(y|x_{0})$ is Gaussian, and the latter is the
score function evaluated at $x_{0}$ and can be replaced with
$S_{\theta}(x_{0},0)$.

In practice, $S_{\theta}(x_{0},0)$ is only accurate when $x_{0}$ lies closer
to the data manifold and is typically inaccurate for $x_{0}$ in low-likelihood
regions outside the data manifold (refer to [23]). This makes the score
estimate inaccurate at the start (when $x_{0}$ is initialized randomly) and
also during the gradient ascent since the intermediate $x^{k}_{0}$ are not
constrained to lie on the data manifold. This issue is avoided when inverse
problems are typically solved through the reverse diffusion process (Sec. 2.2)
which drifts the noisy sample towards the data manifold using
$S_{\theta}(x_{t},t)$ while simultaneously ensuring measurement consistency.
But there the bottleneck is to estimate $\nabla_{x_{t}}\log P(y|x_{t})$ which
is intractable, as discussed earlier.

Here, we present our proposed MAP formulation. Let $z\sim
P(x_{T})=\mathcal{N}(0,\sigma^{2}(T)\mathbb{I})$, denote a purely noisy
sample, and $\mathcal{M}$ denote the data Manifold. The PF ODE trajectory maps
$z$ to a sample $x_{0}\in\mathcal{M}$, given by $x_{0}=f_{\theta}(z,T)$, where
$f_{\theta}$ is the consistency model. It is also evident that, $\forall\
x_{0}\in\mathcal{M}$, $\exists\ z\sim P(x_{T})$ such that
$x_{0}=f_{\theta}(z,T)$. Hence, the usual MAP formulation in Eq. 8 is
equivalent to the proposed MAP formulation in Eqs. 10 and 11.

$z^{*}=\arg\max_{z}\ \log P(x_{0}=f_{\theta}(z,T)|y)$ (10)
$x^{*}_{0}=f_{\theta}(z^{*},T)$ (11)

With our proposed formulation, we update $z$ with gradient ascent steps as in
Eq. 12 for finding $z^{*}$. The update step now requires computing the
gradient term $\nabla_{z}\log P(f_{\theta}(z,T)|y)$, which can be reformulated
as a vector-jacobian product (vjp) as shown in Eq. 13.

The vector in the vjp from Eq. 13, is the gradient term $\nabla_{x_{0}}\log
P(x_{0}|y)$ evaluated at $x_{0}=f_{\theta}(z,T)$ which lies on the data
manifold (by definition of consistency model), and can be accurately evaluated
unlike the previous case.

$z^{k+1}=z^{k}+\lambda*\nabla_{z}\log P(f_{\theta}(z,T)|y)$ (12)

$\nabla_{z}\log P(f_{\theta}(z,T)|y)={\left(\frac{\partial
f_{\theta}(z,T)}{\partial z}\right)}^{\intercal}\nabla_{x_{0}}\log
P(x_{0}|y)\bigg{|}_{x_{0}=f_{\theta}(z,T)}$

(13)

Here we provide the high-level overview of a practical implementation using
our MAP formulation. In practice, even a consistency model $f_{\theta}$ can
benefit from multi-step sampling [22]. Therefore we propose a multi-step
gradient ascent scheme called _MAP-Gradient-Ascent (MAP-GA)_ as described in
Algorithm 1. Note that $\tau$ refers to a time-step schedule with
$T=\tau_{n}>\tau_{n-1}>..>\tau_{1}>\tau_{0}=0$ and $\sigma$ refers to the
monotonically increasing noise schedule for $t\in[0,T]$, with $\sigma_{0}=0,$
and $\sigma_{T}=\infty$ (high value in practice), $y$ is the measurement,
$f_{\theta}$ and $S_{\theta}$ are the consistency model and the score model
respectively. $num\\_iter$ denotes the number of gradient ascent iterations
per time step, and $\lambda$ denotes the learning rate. Algorithm 1 can be
applied to solve any inverse problem effectively, given that we know the
optimal hyper-parameters, such as the time-step schedule $\tau$, the learning
rate $\lambda$, learning rate schedule, $num\\_steps$, etc. However, finding
those in practice can be quite challenging.

input : $\tau=(\tau_{n},..\tau_{1},\tau_{0})$, _$f_{\theta}$_ , _$S_{\theta}$_
, _y_ , _num_iter_ , _$\lambda$_ , _$\sigma$_

$z\sim\mathcal{N}(0,\sigma^{2}_{\tau_{n}}\mathbb{I})$

for  _$i\ in\ (n,n-1,..,1)$_ do

$t=\tau_{i}$

for  _$j\ in\ (1,2,..,num\\_iter)$_ do

$z=z+\lambda*\nabla_{z}\log P(f_{\theta}(z,t)|y)$

end for

$\hat{x}_{0}=f_{\theta}(z,t)$

$z=\mathcal{N}(\hat{x}_{0},\sigma^{2}_{\tau_{i-1}}\mathbb{I})$

end for

output : z

Algorithm 1 MAP-GA (MAP-Gradient-Ascent)

###  3.3 Practical Implementation

In practice, to avoid numerical issues, instead of solving for $x^{*}_{0}$, we
solve for $x^{*}_{\epsilon}=\arg\max_{x_{\epsilon}}\ \log P(x_{\epsilon}|y)$,
for a small $\epsilon$ such that $\sigma_{\epsilon}\approx 0$. So we solve the
MAP formulation in Eqs. 14 and 15. Since,
$P(x_{\epsilon}|x_{0})=\mathcal{N}(x_{0},\sigma^{2}_{\epsilon}\mathbb{I})$,
for very small values of $\sigma_{\epsilon}$, the distinction between $x_{0}$
and $x_{\epsilon}$ remain insignificant for all practical purposes. The
consistency models in [22] are also learned to map the points on PF ODE
trajectory to the corresponding $x_{\epsilon}$ instead of $x_{0}$.

$z^{*}=\arg\max_{z}\ \log P(x_{\epsilon}=f_{\theta}(z,T)|y)$ (14)
$x^{*}_{0}\approx x^{*}_{\epsilon}=f_{\theta}(z^{*},T)$ (15)

Here, we describe in detail the computation of the gradient term
$\nabla_{z}\log P(f_{\theta}(z,t)|y)$ in practice. From Eqs. 16 and 17, this
requires the estimation of the gradient of log-likelihood i.e
$\nabla_{x_{\epsilon}}\log P(y|x_{\epsilon})$ and the gradient of log-prior
i.e $\nabla_{x_{\epsilon}}\log P(x_{\epsilon})$. The terms $P(y|x_{\epsilon})$
and $P(x_{\epsilon})$ are also referred to as the likelihood and the prior
respectively.

$\nabla_{z}\log P(f_{\theta}(z,t)|y)={\left(\frac{\partial
f_{\theta}(z,t)}{\partial z}\right)}^{\intercal}\nabla_{x_{\epsilon}}\log
P(x_{\epsilon}|y)\bigg{|}_{x_{\epsilon}=f_{\theta}(z,t)}$

(16)

$\nabla_{x_{\epsilon}}\log
P(x_{\epsilon}|y)\bigg{|}_{x_{\epsilon}=f_{\theta}(z,t)}=\biggl{\\{}\nabla_{x_{\epsilon}}\log
P(y|x_{\epsilon})+\nabla_{x_{\epsilon}}\log
P(x_{\epsilon})\biggr{\\}}\bigg{|}_{x_{\epsilon}=f_{\theta}(z,t)}$

(17)

Computing the gradient of log-likelihood
$P(y|x_{0})=\mathcal{N}(\mathcal{A}(x_{0}),\sigma^{2}_{y}\mathbb{I})$, and
$P(x_{\epsilon}|x_{0})=\mathcal{N}(x_{0},\sigma^{2}_{\epsilon}\mathbb{I})$,
given by Eq. 1 and the diffusion perturbation kernel respectively. Since
$\sigma_{\epsilon}\approx 0$, $P(x_{0}|x_{\epsilon})$ can be approximated with
$P(x_{\epsilon}|x_{0})$, and for a linear forward operator i.e
$\mathcal{A}=H\in\mathbb{R}^{m\times n}$, we can derive
$P(y|x_{\epsilon})=\mathcal{N}(Hx_{\epsilon},\sigma^{2}_{y}\mathbb{I}+\sigma^{2}_{\epsilon}HH^{\intercal})$.
For non-linear $\mathcal{A}$, by linearizing it around $x_{\epsilon}$, similar
approximations can be derived. As $\sigma_{\epsilon}\to 0$, these
approximations tend to hold well. Given the tractable form of the likelihood
term above, computing the gradient of the log-likelihood is apparent.

Computing the gradient of log-prior
The gradient of the log prior i.e $\nabla_{x_{\epsilon}}\log P(x_{\epsilon})$
is essentially the score function evaluated at $x_{\epsilon}$. Given a score
function $S_{\theta}(x,t)$, learned by the unconditional diffusion model,
$\nabla_{x_{\epsilon}}\log P(x_{\epsilon})=S_{\theta}(x_{\epsilon},\epsilon)$.
Learning the score function is equivalent to learning a denoiser, and vice-
versa [10]. Hereon, we denote the unconditional diffusion model as learning
the denoiser $D_{\theta}(x,t)$, from which the score function can be computed
as  $\nabla_{x_{t}}\log P(x_{t})$
$=\frac{D_{\theta}(x_{t},t)-x_{t}}{\sigma^{2}_{t}}$.

##  4 Image Inpainting with MAP-GA

In Computer vision, image inpainting [30, 14] is a typical image restoration
task. Given a masked image (and the corresponding binary mask), the goal here
is to recover (reconstruct) the missing pixels of the masked image. Image
inpainting can be modeled as a linear inverse problem. Let
$x\in\mathbb{R}^{n}$ denote the original image, $y\in\mathbb{R}^{m}$ denote a
masked image with only visible pixels ($m\leq n$) and $H\in\mathbb{R}^{m\times
n}$ denotes a corresponding forward operator for a given mask. The inpainting
problem is characterized by $y=Hx+\eta$, where
$\eta\sim\mathcal{N}(0,\sigma^{2}_{y}\mathbb{I})$. Note that a given mask
defines $H$ and for inpainting, $H$ has a defined structure, where the rows of
$H$ are one-hot and are orthogonal. This also implies
$HH^{\intercal}=\mathbb{I}_{m\times m}$.

input : time schedule: $\tau=[\ \tau_{n},\tau_{n-1},..\tau_{1},\tau_{0}\ ]$,

noise schedule: $\sigma(.)$, Denoiser: $D_{\theta},$

consistency model: $C_{\theta}$,

measurement: $y$, learning rate: $\lambda,$

num gradient ascent iter: num_iter,

boolean: use_prior (default True),

Inpainting matrix: H,

Measurement noise: $\sigma_{y}$

(Note: $HH^{\intercal}=\mathbb{I}_{m\times m}$, $\tau_{0}=\epsilon$,
$\tau_{n}=T$),

$z\sim\mathcal{N}(0,\sigma^{2}_{\tau_{n}}\mathbb{I})$

for  _$i\ in\ (\ n,n-1,\ ..\ 1\ )$_ do

$t=\tau_{i}$

for  _$j\ in\ (\ 1,2,\ ..\ num\\_iter\ )$_ do

$\hat{x}_{\epsilon}=C_{\theta}(z,t)$

$grad_{likelihood}=\frac{H^{\intercal}(y-H\hat{x}_{\epsilon})}{\sigma^{2}_{y}+\sigma^{2}_{\epsilon}}$

if  _use_prior_ then

$grad_{prior}=\frac{D_{\theta}(\hat{x}_{\epsilon},\epsilon)-\hat{x}_{\epsilon}}{\sigma^{2}_{\epsilon}}$

end if

else

$grad_{prior}=0$

end if

$grad_{posterior}=grad_{likelihood}+grad_{prior}$

$grad=\left(\frac{\partial C_{\theta}(z,t)}{\partial
z}\right)^{\intercal}grad_{posterior}$

$z=z+\lambda*grad$

end for

$\hat{x}_{\epsilon}=C_{\theta}(z,t)$

$z=\mathcal{N}(\hat{x}_{\epsilon},\sigma^{2}_{\tau_{i-1}}-\sigma^{2}_{\tau_{0}}\mathbb{I})$

end for

output : z

Algorithm 2 MAP-GA (Algorithm 1 expanded for Inpainting)

In this work, to show the effectiveness of our proposed MAP framework in
practice, we consider the task of image inpainting with diverse mask settings.
We expand Algorithm 1 for the task of image inpainting, and include all the
specific details in Algorithm 2. The core quantity in the algorithm that needs
to be evaluated is $\nabla_{z}\log P(f_{\theta}(z,t)|y)$, which involves the
estimation of gradients of the log-likelihood and the log-prior terms (Eqs. 16
and 17). In the algorithm, we make a choice (indicated by the $use\\_prior$
keyword) of retaining or dropping the gradient of the log-prior. Dropping the
prior term implies a choice of uniform prior, and the algorithm now optimizes
for the maximum likelihood estimate instead of the MAP estimate, by
considering any sample (which is consistent with our measurement) on the data
manifold to be equally good. This is also the setting considered in the
concurrent works [4, 28].

Note that the algorithm makes use of both the consistency model ($C_{\theta}$)
and the denoiser ($D_{\theta}$), but this makes the method more demanding
compared to other methods that only use the denoiser. For that, we make an
argument as follows. We use the pre-trained denoiser and the consistency model
from [10] and [22] respectively, with noise schedule $\sigma(t)=t$ and
$t\in[\epsilon,T]$, where $\epsilon=0.002,T=80$. Consider the PF ODE for the
above setting, as follows.

$\displaystyle dx$ $\displaystyle=-\sigma(t)\dot{\sigma}(t)\nabla_{x_{t}}\log
P(x_{t})dt$
$\displaystyle=-\sigma(t)\dot{\sigma}(t)\frac{D_{\theta}(x,t)-x}{\sigma^{2}(t)}dt$
$\displaystyle=-\dot{\sigma}(t)\frac{D_{\theta}(x,t)-x}{\sigma(t)}dt$
$\displaystyle=-\frac{D_{\theta}(x,t)-x}{t}dt$

This PF ODE determines the trajectory and solving the trajectory origin
$x_{0}$ involves solving the ODE above. Note that this $x_{0}$ is what the
consistency model $C_{\theta}$ is trained to predict. As a rough
approximation, solving the ODE with a backward Euler discretization step from
$t$ to $0$, which essentially assumes the trajectory curves are linear, i.e.
the jacobian $\frac{dx}{dt}$ is constant for the interval $[0,t]$, and this
gives $C_{\theta}(x_{t},t)=x_{0}\approx D_{\theta}(x_{t},t)$, with the
approximation getting more accurate as the trajectory curve gets more linear.
In an empirical setting, [10] also observes trajectories become more linear
when $\sigma(t)\to 0$. While more analysis on this is still due, it motivates
us to look at the denoiser as a proxy of the consistency model, which
gradually becomes more and more accurate as $\sigma(t)\to 0$. Hence, we also
consider settings that replace the consistency model with the denoiser in our
algorithm.

Figure 1: From left to right: original image, and mask settings: box50, half,
expand, box25, sr2x, altlines

input : noise schedule: $\sigma(.)$, Denoiser: $D_{\theta},$

consistency model: $C_{\theta}$,

measurement: $y$, learning rate: $\lambda,$

num gradient ascent iter: num_iter,

Inpainting matrix: H,

Measurement noise: $\sigma_{y}$

$\tau^{map-
ga}=[\tau^{map}_{n},\tau^{map}_{n-1},..\tau^{map}_{1},\tau^{map}_{0}]$,

$\tau^{pgdm}=[\tau^{pgdm}_{m},\tau^{pgdm}_{m-1},..\tau^{pgdm}_{1},\tau^{pgdm}_{0}]$,

(Note: $HH^{\intercal}=\mathbb{I}_{m\times m}$)

(Note: $\sigma_{\tau^{map}_{0}}=\sigma_{y}$ and $\tau^{map}_{n}=T$)

(Note: $\sigma_{\tau^{pgdm}_{m}}=\sigma_{y}$ and $\tau^{pgdm}_{0}=\epsilon$),

=======MAP-GA=======

$z\sim\mathcal{N}(0,\sigma^{2}_{\tau^{map}_{n}}\mathbb{I})$

for  _$i\ in\ (\ n,n-1,\ ..\ 1\ )$_ do

$t=\tau^{map}_{i}$

for  _$j\ in\ (\ 1,2,\ ..\ num\\_iter\ )$_ do

$\hat{x}_{\epsilon}=C_{\theta}(z,t)$

$grad_{likelihood}=\frac{H^{\intercal}(y-H\hat{x}_{\epsilon})}{\sigma^{2}_{y}+\sigma^{2}_{\epsilon}}$

$grad_{posterior}=grad_{likelihood}$

$grad=\left(\frac{\partial C_{\theta}(z,t)}{\partial
z}\right)^{\intercal}grad_{posterior}$

$z=z+\lambda*grad$

end for

$\hat{x}_{\epsilon}=C_{\theta}(z,t)$

$z=\mathcal{N}(\hat{x}_{\epsilon},\sigma^{2}_{\tau^{map}_{i-1}}-\sigma^{2}_{\tau^{map}_{0}}\mathbb{I})$

end for

========PGDM========

$x_{\tau^{pgdm}_{m}}=z$

for  _$i\ in\ (\ m,m-1,\ ..\ 1\ )$_ do

$t=\tau^{pgdm}_{i}$

$\hat{x}_{t}=D_{\theta}(z,t)$

$\mu_{t}=\hat{x}_{t}+\sigma^{2}_{t}*\nabla_{x_{t}}\log P(y|x_{t})$ *pgdm
update*

$x_{\tau^{pgdm}_{i-1}}=\mathcal{N}(\mu_{t},\sigma^{2}_{\tau^{pgdm}_{i-1}}-\sigma^{2}_{\tau^{pgdm}_{0}}\mathbb{I})$

end for

output : $x_{\tau^{pgdm}_{0}}$

Algorithm 3 MAP-GA-PGDM for Inpainting

When using Algorithm 2 for noisy image inpainting ($\sigma_{y}>0$) in
practice, we observe that it requires careful tuning of the learning rate and
other hyperparameters. To avoid such additional and sensitive hyperparameters,
we present Algorithm 3 for noisy inpainting, based on our empirical
observations in Tab. 2. The algorithm requires no additional hyperparameters
and works well in practice. The idea behind Algorithm 3, is to find a good
approximate solution using MAP-GA at diffusion time $t=\tau$ where
$\sigma_{\tau}=\sigma_{y}$, and then use it as an initialization for PGDM [20]
for $\sigma_{t}<\sigma_{y}$. More specifically, we use MAP-GA until
$\sigma_{t}=\sigma_{y}$, to estimate $x_{\sigma_{y}}$ and later use this as
initialization to PGDM for $\sigma_{t}<\sigma_{y}$. (Note that we could use
$t$ and $\sigma_{t}$ interchangeably in our notations since our noise schedule
is $\sigma(t)=t$). We do not use the prior term for the MAP-GA part in
Algorithm 3, as we find it more effective.

##  5 Experiments

In all our experiments, we use the pre-trained class-conditional denoiser from
[10] and the consistency model from [22] with their default settings. Inspired
by [12], we evaluate the inpainting performance across six different mask
settings that denote varying levels of degradation. Examples of our mask
settings can be seen in Fig. 1. The mask settings box50 and box25 indicate a
square crop at the center of the image, with the crop width equal to 50% and
25% of the image width respectively. In half, we mask out the right half of
the image, the mask expand is the complement of box25. sr2x denotes a 2x
super-resolution mask, and altlines masks out alternate rows of pixels. We
evaluate the inpainted images using FID and LPIPS metrics.

In Tab. 1, we evaluate the performance of MAP-GA (Algorithm 2) for noiseless
image inpainting (i.e. $\sigma_{y}=0$) on Imagenet [17] validation set with
64x64 resolution having 1k images. In the experiments, the setting MAP-GA
denotes the default Algorithm 2, MAP-GA(NP) denotes MAP-GA with no prior, MAP-
GA(D) denotes MAP-GA with denoiser replacing the consistency model, and MAP-
GA(D,NP) denotes MAP-GA with no prior and the denoiser replacing the
consistency model. The results from Tab. 1 show MAP-GA and variants
outperforming PGDM by a significant margin. Note that MAP-GA uses gradient
ascent (first-order gradient-based optimization algorithm) to optimize the
underlying MAP objective which is typically highly non-convex and is not
guaranteed to find the global optima. However, with an initialization, closer
to the global optima, MAP-GA could converge.

| box50 | half | expand | box25 | sr2x | altlines
---|---|---|---|---|---|---
Method |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$
|  |  |  |  |  |  |  |  |  |  |  |
MAP-GA(D,NP) | 32.312 | 0.096 | 37.986 | 0.157 | $93.068$ | $0.419$ | 11.019 | 0.019 | 30.072 | 0.034 | 18.282 | 0.015
MAP-GA(D) | $35.761$ | $0.100$ | $43.315$ | $0.164$ | $106.32$ | $0.430$ | $11.169$ | $0.019$ | $33.188$ | $0.036$ | $19.103$ | $0.016$
MAP-GA(NP) | $34.944$ | $0.113$ | $39.243$ | $0.162$ | 69.004 | 0.388 | $12.818$ | $0.027$ | $30.303$ | $0.035$ | $20.321$ | $0.018$
MAP-GA | $36.733$ | $0.113$ | $41.151$ | $0.164$ | $87.952$ | $0.400$ | $12.836$ | $0.027$ | $33.752$ | $0.036$ | $20.895$ | $0.018$
PGDM [20] | $49.370$ | $0.151$ | $54.261$ | $0.245$ | $127.95$ | $0.479$ | $14.255$ | $0.021$ | $38.433$ | $0.046$ | $20.446$ | $0.019$

Table 1:  Noiseless inpainting on ImageNet64 1k validation set using MAP-GA
and variants. The setting MAP-GA denotes the default Algorithm 2, MAP-GA(NP)
denotes MAP-GA with no prior, MAP-GA(D) denotes MAP-GA with denoiser replacing
the consistency model, and MAP-GA(D,NP) denotes MAP-GA with no prior and the
denoiser replacing the consistency model.

| box50 | half | expand | box25 | sr2x | altlines
---|---|---|---|---|---|---
Method |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$
|  |  |  |  |  |  |  |  |  |  |  |
MAP-GA(D,NP) | $26.661$ | $0.043$ | $30.621$ | $0.064$ | $73.510$ | $0.155$ | 9.165 | $0.010$ | $29.878$ | $0.027$ | 16.971 | 0.012
MAP-GA(D) | $34.689$ | $0.052$ | $41.098$ | $0.083$ | $90.768$ | $0.175$ | $10.891$ | $0.011$ | $34.334$ | $0.031$ | $18.718$ | $0.013$
MAP-GA(NP) | 25.205 | 0.043 | 27.208 | 0.047 | 42.671 | 0.075 | $10.936$ | $0.017$ | 27.496 | 0.023 | $18.732$ | $0.013$
MAP-GA | $29.508$ | $0.045$ | $33.653$ | $0.050$ | $73.318$ | $0.106$ | $11.257$ | $0.016$ | $34.364$ | $0.026$ | $20.041$ | $0.013$
PGDM [20] | $31.952$ | $0.056$ | $36.639$ | $0.082$ | $71.352$ | $0.138$ | $10.886$ | 0.009 | $33.199$ | $0.041$ | $19.551$ | $0.018$

Table 2:  Noiseless inpainting on ImageNet64 1k validation set. Using the
ground truth image ($x_{0}$) for the measurement $y$, we create a sample at
$t=0.5$ via ($x_{0.5}=x_{0}+0.5*\eta$, where,
$\eta\sim\mathcal{N}(0,\mathbb{I}))$ and initialize Algorithm 2 with
$z=x_{0.5}$, and $\tau_{n}=0.5$

To corroborate this, we design the following toy experiment. We consider the
same noiseless inpainting task from earlier, but, in Algorithm 2, we set
$\tau_{n}=\hat{\tau}\ll T$, instead of the default setting $\tau_{n}=T$ and we
also initialize $z\sim\mathcal{N}(x_{0},\sigma^{2}_{\hat{\tau}}\mathbb{I})$,
where $x_{0}$ is the corresponding ground truth image for the measurement $y$.
From Tab. 2, we observe significant improvements in the performance. (We use
$\hat{\tau}=0.5$, note that $\sigma_{0.5}=0.5$, as we use the noise schedule
$\sigma(t)=t$). This shows that MAP-GA is only limited by the choice of
optimizer and reinforces the need for better optimization algorithms. While it
is important to consider better design choices (for the schedules and
hyperparameters), adaptive-gradient-based optimizers (such as momentum, adam),
or higher-order methods, we leave this for future work as it requires a
thorough analysis.

In Tab. 3, we report the results on noisy inpainting using Algorithm 3 on
Imagenet 1k validation set with 64x64 resolution. $\sigma_{y}$ indicates the
measurement noise and the setting MAP-GA-PGDM denote the default Algorithm 3,
MAP-GA-PGDM(D) denote MAP-GA-PGDM with denoiser replacing the consistency
model. In Algorithm 3, and all its variants, we fix $m=20$ time steps for
PGDM. Even with a high level of measurement noise, MAP-GA-PGDM shows promising
improvements over PGDM.

In all our experiments, each of our methods has been given a budget of 1000
steps, and we run ablations for (num_time_steps, num_iter) from the set
{(20,50),(50,20),(100,10),(200,5),(250,4),(500,2),(1000,1)}. We follow a time-
step schedule similar to [10], with their default settings. The best results
for all our methods resulted from the ablation (20,50) denoting 20 time-steps
with 50 gradient ascent iterations. For PGDM [20], each setting from our
results table is run with different choices of num_steps from
{20,50,100,200,250,500,1000}, and most of the runs have their peak performance
at $500$ steps. The performance reported in the tables is the best result for
each method over the ablations for (num_time_steps, num_iter). The choice of
learning rate in all our experiments was
$\sigma^{2}_{y}+\sigma^{2}_{\epsilon}$.

|  | box50 | half | expand | box25 | sr2x | altlines
---|---|---|---|---|---|---|---
Method | $\sigma_{y}$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$ |  FID$\downarrow$ |  LPIPS$\downarrow$
|  |  |  |  |  |  |  |  |  |  |  |  |
MAP-GA-PGDM(D) | $0.05$ | 56.173 | 0.125 | $62.026$ | $0.193$ | $108.16$ | $0.438$ | 38.725 | 0.039 | $57.41$ | $0.080$ | $46.046$ | $0.046$
MAP-GA-PGDM | $0.05$ | $58.283$ | $0.147$ | 61.588 | 0.184 | 91.508 | 0.406 | $44.667$ | $0.085$ | 57.308 | 0.073 | 45.265 | 0.042
PGDM [20] | $0.05$ | $77.824$ | $0.175$ | $80.248$ | $0.257$ | $135.99$ | $0.495$ | $53.136$ | $0.049$ | $86.289$ | $0.126$ | $66.203$ | $0.066$
|  |  |  |  |  |  |  |  |  |  |  |  |
MAP-GA-PGDM(D) | $0.1$ | 72.543 | 0.166 | $79.535$ | $0.230$ | $114.51$ | $0.464$ | 57.556 | 0.076 | $76.733$ | $0.145$ | $65.154$ | $0.096$
MAP-GA-PGDM | $0.1$ | $74.130$ | $0.191$ | 78.322 | 0.225 | 103.85 | 0.440 | $65.720$ | $0.161$ | 76.248 | 0.134 | 63.134 | 0.089
PGDM [20] | $0.1$ | $96.485$ | $0.216$ | $99.170$ | $0.286$ | $145.40$ | $0.519$ | $78.620$ | $0.100$ | $109.47$ | $0.231$ | $90.925$ | $0.138$
|  |  |  |  |  |  |  |  |  |  |  |  |
MAP-GA-PGDM(D) | $0.2$ | $91.112$ | 0.243 | $99.432$ | $0.297$ | $119.45$ | $0.496$ | 79.329 | 0.154 | $110.48$ | $0.260$ | $87.727$ | $0.188$
MAP-GA-PGDM | $0.2$ | 91.087 | $0.252$ | 97.531 | 0.296 | 116.78 | 0.491 | $80.647$ | $0.177$ | 103.66 | 0.248 | 86.174 | 0.186
PGDM [20] | $0.2$ | $116.81$ | $0.292$ | $117.03$ | $0.344$ | $159.17$ | $0.562$ | $103.43$ | $0.193$ | $132.89$ | $0.366$ | $114.54$ | $0.255$

Table 3:  Noisy inpainting on ImageNet64 1k validation set. $\sigma_{y}$
denotes the measurement noise. The setting MAP-GA-PGDM denotes the default
Algorithm 3, MAP-GA-PGDM(D) denote MAP-GA-PGDM with denoiser replacing the
consistency model.

##  6 Discussion

The algorithms presented in this paper are similar to the concurrent works
ZSIR [4], and DMPlug [28] from a practical perspective. However, our proposed
MAP formulation has a strong theoretical motivation and connects the
probability flow ODE and the consistency model with the MAP optimization for
solving inverse problems. While we compare only against PGDM [20], it is the
most competitive baseline and can also be seen from the experiments in DMPlug
[28]. In this work, we analyze our method extensively, on image inpainting
with various degradation masks. However, the experiments from DMPlug [28]
further complement and validate our work. We show that the gradient of the
log-prior is tractable, making the gradient of log-posterior tractable.
However, [4, 28] considers an arbitrary choice of uniform prior and ignores
the prior term altogether. We show that MAP-GA is only limited by the choice
of the optimizer in practice and our proposed framework has several directions
for future research, such as improving optimization with hessian-based
methods, more analysis on local optima, theory connecting the denoiser with
the consistency model among others. Also, our method and algorithms are valid
for general SDE-based Diffusion models and for more general inverse problems
than the inpainting task we considered.

Figure 2: Noiseless inpainting example. From left to right: original image,
masked image, MAP-GA(Algorithm 2), PGDM($\sigma_{y}=0$). From top to bottom:
_half mask_ , _expand mask_ , _box50 mask_

Figure 3: Noisy inpainting example. From left to right: original image, masked
image, MAP-GA-PGDM(Algorithm 3, $\sigma_{y}=0.1$), PGDM($\sigma_{y}=0.1$).
From top to bottom: _sr2x mask_ , _box25 mask_

##  7 Conclusion

In this paper, we proposed a novel MAP formulation for solving inverse
problems using pre-trained unconditional diffusion models. Note that
conditional generation is a core requirement in solving inverse problems. We
connect the probability flow ODE and the consistency model with the
optimization process for the MAP objective in tasks that involve conditional
generation. We show that the gradient of the MAP objective is tractable,
allowing the use of gradient-based optimization methods. To use our framework
in practice, we propose an algorithm with a multi-step gradient ascent
strategy for MAP optimization. Since we use gradient-ascent (a simple first-
order gradient-based method), we observed that it requires extensive tuning of
the learning rate and other hyperparameters due to the challenging
optimization objective. To address this, we propose practically effective
algorithms, for image inpainting using our framework (which can be extended to
other inverse problems). We validate our framework with extensive experiments
on image inpainting across diverse mask settings.

## Acknowledgements

This work is part of the Marie Skłodowska-Curie Actions project _MODELAIR_ ,
funded by the European Commission under the Horizon Europe program through
grant agreement no. 101072559. The computations and the data handling were
enabled by resources provided by the National Academic Infrastructure for
Supercomputing in Sweden (NAISS), partially funded by the Swedish Research
Council through grant agreement no. 2022-06725. Bharath thanks Sebastian
Gerard and Heng Fang for their feedback on improving the presentation of this
paper.

## References

  * [1] Simon Arridge, Peter Maass, Ozan Öktem, and Carola-Bibiane Schönlieb.  Solving inverse problems using data-driven models.  Acta Numerica, 28:1–174, 2019.
  * [2] Ashish Bora, Ajil Jalal, Eric Price, and Alexandros G Dimakis.  Compressed sensing using generative models.  In International conference on machine learning, pages 537–546. PMLR, 2017.
  * [3] Benjamin Boys, Mark Girolami, Jakiw Pidstrigach, Sebastian Reich, Alan Mosca, and O Deniz Akyildiz.  Tweedie moment projected diffusions for inverse problems.  arXiv preprint arXiv:2310.06721, 2023.
  * [4] Hamadi Chihaoui, Abdelhak Lemkhenter, and Paolo Favaro.  Zero-shot image restoration via diffusion inversion, 2024.
  * [5] Hyungjin Chung, Jeongsol Kim, Michael Thompson Mccann, Marc Louis Klasky, and Jong Chul Ye.  Diffusion posterior sampling for general noisy inverse problems.  In International Conference on Learning Representations, 2023.
  * [6] Hyungjin Chung, Suhyeon Lee, and Jong Chul Ye.  Decomposed diffusion sampler for accelerating large-scale inverse problems.  arXiv preprint arXiv:2303.05754, 2023.
  * [7] Hyungjin Chung, Byeongsu Sim, Dohoon Ryu, and Jong Chul Ye.  Improving diffusion models for inverse problems using manifold constraints.  Advances in Neural Information Processing Systems, 35:25683–25696, 2022.
  * [8] Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio.  Generative adversarial nets.  Advances in neural information processing systems, 27, 2014.
  * [9] Jonathan Ho, Ajay Jain, and Pieter Abbeel.  Denoising diffusion probabilistic models.  Advances in neural information processing systems, 33:6840–6851, 2020.
  * [10] Tero Karras, Miika Aittala, Timo Aila, and Samuli Laine.  Elucidating the design space of diffusion-based generative models.  Advances in neural information processing systems, 35:26565–26577, 2022.
  * [11] Bahjat Kawar, Michael Elad, Stefano Ermon, and Jiaming Song.  Denoising diffusion restoration models.  Advances in Neural Information Processing Systems, 35:23593–23606, 2022.
  * [12] Andreas Lugmayr, Martin Danelljan, Andres Romero, Fisher Yu, Radu Timofte, and Luc Van Gool.  Repaint: Inpainting using denoising diffusion probabilistic models.  In Proceedings of the IEEE/CVF conference on computer vision and pattern recognition, pages 11461–11471, 2022.
  * [13] Alex Nichol, Prafulla Dhariwal, Aditya Ramesh, Pranav Shyam, Pamela Mishkin, Bob McGrew, Ilya Sutskever, and Mark Chen.  Glide: Towards photorealistic image generation and editing with text-guided diffusion models.  arXiv preprint arXiv:2112.10741, 2021.
  * [14] Deepak Pathak, Philipp Krahenbuhl, Jeff Donahue, Trevor Darrell, and Alexei A Efros.  Context encoders: Feature learning by inpainting.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2536–2544, 2016.
  * [15] Xinyu Peng, Ziyang Zheng, Wenrui Dai, Nuoqian Xiao, Chenglin Li, Junni Zou, and Hongkai Xiong.  Improving diffusion models for inverse problems using optimal posterior covariance.  In Forty-first International Conference on Machine Learning, 2024.
  * [16] Danilo Rezende and Shakir Mohamed.  Variational inference with normalizing flows.  In International conference on machine learning, pages 1530–1538. PMLR, 2015.
  * [17] Olga Russakovsky, Jia Deng, Hao Su, Jonathan Krause, Sanjeev Satheesh, Sean Ma, Zhiheng Huang, Andrej Karpathy, Aditya Khosla, Michael Bernstein, et al.  Imagenet large scale visual recognition challenge.  International journal of computer vision, 115:211–252, 2015.
  * [18] Chitwan Saharia, William Chan, Huiwen Chang, Chris Lee, Jonathan Ho, Tim Salimans, David Fleet, and Mohammad Norouzi.  Palette: Image-to-image diffusion models.  In ACM SIGGRAPH 2022 conference proceedings, pages 1–10, 2022.
  * [19] Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli.  Deep unsupervised learning using nonequilibrium thermodynamics.  In International conference on machine learning, pages 2256–2265. PMLR, 2015.
  * [20] Jiaming Song, Arash Vahdat, Morteza Mardani, and Jan Kautz.  Pseudoinverse-guided diffusion models for inverse problems.  In International Conference on Learning Representations, 2023.
  * [21] Yang Song and Prafulla Dhariwal.  Improved techniques for training consistency models.  arXiv preprint arXiv:2310.14189, 2023.
  * [22] Yang Song, Prafulla Dhariwal, Mark Chen, and Ilya Sutskever.  Consistency models.  arXiv preprint arXiv:2303.01469, 2023.
  * [23] Yang Song and Stefano Ermon.  Generative modeling by estimating gradients of the data distribution.  Advances in neural information processing systems, 32, 2019.
  * [24] Yang Song, Sahaj Garg, Jiaxin Shi, and Stefano Ermon.  Sliced score matching: A scalable approach to density and score estimation.  In Uncertainty in Artificial Intelligence, pages 574–584. PMLR, 2020.
  * [25] Yang Song, Jascha Sohl-Dickstein, Diederik P Kingma, Abhishek Kumar, Stefano Ermon, and Ben Poole.  Score-based generative modeling through stochastic differential equations.  arXiv preprint arXiv:2011.13456, 2020.
  * [26] Pascal Vincent.  A connection between score matching and denoising autoencoders.  Neural computation, 23(7):1661–1674, 2011.
  * [27] Ricardo Vinuesa and Steven L Brunton.  Enhancing computational fluid dynamics with machine learning.  Nature Computational Science, 2(6):358–366, 2022.
  * [28] Hengkang Wang, Xu Zhang, Taihui Li, Yuxiang Wan, Tiancong Chen, and Ju Sun.  Dmplug: A plug-in method for solving inverse problems with diffusion models.  arXiv preprint arXiv:2405.16749, 2024.
  * [29] Yinhuai Wang, Jiwen Yu, and Jian Zhang.  Zero-shot image restoration using denoising diffusion null-space model.  arXiv preprint arXiv:2212.00490, 2022.
  * [30] Raymond A Yeh, Chen Chen, Teck Yian Lim, Alexander G Schwing, Mark Hasegawa-Johnson, and Minh N Do.  Semantic image inpainting with deep generative models.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 5485–5493, 2017.
  * [31] Yuanzhi Zhu, Kai Zhang, Jingyun Liang, Jiezhang Cao, Bihan Wen, Radu Timofte, and Luc Van Gool.  Denoising diffusion models for plug-and-play image restoration.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1219–1229, 2023.