# GoMVS: Geometrically Consistent Cost Aggregation for Multi-View Stereo

Jiang Wu1 * Rui Li1,2 * Haofei Xu2,3 Wenxun Zhao1 Yu Zhu1 ✝ Jinqiu Sun1
Yanning Zhang1✝
1Northwestern Polytechnical University 2ETH Zürich 3University of Tübingen,
Tübingen AI Center

###### Abstract

Matching cost aggregation plays a fundamental role in learning-based multi-
view stereo networks. However, directly aggregating adjacent costs can lead to
suboptimal results due to local geometric inconsistency. Related methods
either seek selective aggregation or improve aggregated depth in the 2D space,
both are unable to handle geometric inconsistency in the cost volume
effectively. In this paper, we propose GoMVS to aggregate geometrically
consistent costs, yielding better utilization of adjacent geometries. More
specifically, we correspond and propagate adjacent costs to the reference
pixel by leveraging the local geometric smoothness in conjunction with surface
normals. We achieve this by the geometric consistent propagation (GCP) module.
It computes the correspondence from the adjacent depth hypothesis space to the
reference depth space using surface normals, then uses the correspondence to
propagate adjacent costs to the reference geometry, followed by a convolution
for aggregation. Our method achieves new state-of-the-art performance on DTU,
Tanks & Temple, and ETH3D datasets. Notably, our method ranks 1st on the Tanks
& Temple Advanced benchmark. Code is available at
https://github.com/Wuuu3511/GoMVS.

(a) MVSFormer [2] (b) RA-MVSNet [45] (c) ET-MVSNet [14] (d) GeoMVSNet [46] (e)
Ours

Figure 1: Comparison of reconstruction errors on Tanks and Temple benchmark.
We show precision and recall error maps for the “Horse” scan. Our method
demonstrates notable improvements over existing methods in challenging areas.

†† * indicates equal contributions and ✝ indicates corresponding authors.

##  1 Introduction

Multi-view stereo (MVS) is a fundamental computer vision problem that recovers
3D shapes from posed images by multi-view correspondence matching [21]. Recent
learning-based MVS [38, 25, 11, 30] estimates scene depth from the cost volume
computed by geometric matching, which delivers latent geometric cues crucial
for the final depth [7]. However, the initial cost volume can suffer from
challenging matching conditions, _e.g_., varying illumination, textless areas,
or repetitive patterns, leading to suboptimal pixel-wise costs that hamper
accurate estimations.

To mitigate this issue, cost aggregation plays an important role in removing
matching ambiguities and improving discriminativeness by using the neighboring
information. However, the adjacent costs may deliver inconsistent depth cues
due to the gradual changes in local geometry. As a result, the aggregated
costs are not geometrically guaranteed to have the highest matching score at
the real reference depth, leading to suboptimal depth predictions. The widely
adopted cascade framework [7] can potentially exacerbate this issue as the
adjacent costs can have more divergent costs due to the shifted depth
hypotheses.

As the geometric inconsistency is a common challenge in multi-view stereo and
2-view stereo matching, related methods either adopt learned aggregation [25,
31] or enforce consistency to the aggregated depth [9, 15, 42]. Specifically,
some methods [25, 31] adopt adaptive aggregation schemes to allow networks to
select pixels that potentially correlate well and contribute to the reference
pixel’s geometry. However, they heavily rely on network capabilities and do
not guarantee geometric plausibility from the selected costs. Other methods
[9, 19] seek to refine or regularize the aggregated depth values using jointly
estimated surface normals. However, these methods only refine the output depth
in 2D image space and are inherently unable to handle inconsistencies in the
cost volume, which is vital for MVS methods.

In this paper, we propose GoMVS that aggregates geometrically consistent
costs, allowing better utilization of adjacent geometries. Considering that
the local geometry is usually smooth and exhibits gradual changes, we leverage
the local smoothness to correspond and propagate adjacent costs to the
reference cost. We achieve this by the geometrically consistent aggregation
scheme, which operates on the local convolution window and propagates adjacent
costs with the geometrically consistent propagation (GCP) module. The GCP
module computes the correspondences from the adjacent cost’s hypothesized
depth space to the reference cost’s depth space, using back-projected depth
hypotheses and the surface normal. Then, it propagates the adjacent costs to
the reference by interpolating cost scores with respect to the correspondence.
After propagating adjacent costs within a local window, we aggregate them
using standard convolutions. Unlike previous methods [19, 9, 15] that refine
the predicted depth in the 2D space, our method incorporates geometric
consistency in the cost space, yielding a better utilization of adjacent
geometries. As surface normal is crucial for corresponding and propagating
local costs, we further investigate different choices of normal predictions.
We find appropriately applying off-the-shelf monocular normal models enables
smooth and robust aggregation across datasets. We conduct extensive
experiments to evaluate our method’s effectiveness, and our method achieves
new state-of-the-art on DTU, Tank & Temple, as well as the ETH3D dataset. Our
contributions are summarized as follows:

  * •

We propose GoMVS to aggregate geometrically consistent costs, allowing better
utilization of adjacent geometries.

  * •

We propose a geometrically consistent propagation (GCP) module that allows
geometrically plausible correspondence and propagation in cost space.

  * •

We investigate different choices of normal computation and find that properly
applying the monocular surface normal model performs well across datasets.

##  2 Related Works

###  2.1 Learning-based MVS Methods

Multi-View Stereo (MVS) aims to reconstruct 3D scenes from multiple posed
images. In recent years, learning-based methods have exhibited promising
results. MVSNet [38] uses differentiable homography to construct the cost
volume and employs a 3D U-Net for regularization. Subsequent works improve
this framework in several ways. RNN-based methods [39, 29, 35]and coarse-to-
fine approaches [7, 3, 36, 17, 25] reduce memory consumption through by
designing efficient structures. Another group of methods [4, 12, 14] devises
local or global attention modules to enhance input feature representations.
MVSFormer [2] incorporates an additional pre-trained transformer network,
enhancing the performance of MVS with a powerful feature extractor. However,
it lacks further exploration in terms of geometry. GeoMVSNet [46] utilizes the
coarse depth map to extract additional geometric features. In addition, [34,
44, 27] have designed pixel-wise visibility modules to handle occlusions.

###  2.2 Cost Volume Aggregation

As cost volume is vital for multi-view depth estimation, recent works
introduce different cost aggregation methods to the depth network. NP-CVP-
MVSNet [37] introduces sparse convolution to aggregate matching costs at the
same depth range. WT-MVSNet [12] employs a cost transformer to generate a more
complete and smoother probability volume. GeoMVSNet [46] incorporates the
coarse probability volume to enhance the matching discriminative ability.
While these methods improve the capability of regularization networks, the
local geometric inconsistency of the cost volume still remains and poses
challenges for the final aggregation results.

###  2.3 Normal Assisted Depth Estimation

Surface normal provides rich geometric details and has been widely applied in
recent years to depth estimation tasks. Traditional MVS methods [32, 33, 20]
optimize depth and normal hypotheses simultaneously by constructing a planar
prior model. Inspired by traditional methods, SP-Net [28] performs slanted
plane cost aggregation by learning parameterized local planes. NAPV-MVS [24]
uses local normal similarity to emphasize the most relevant adjacent costs.
NR-MVSNet [10] utilizes depth-normal consistency to adaptively expand the
hypothesis range, providing broader matchings to assist depth inference.
However, these methods do not address the local inconsistency issue. GeoNet
[19] proposes a monocular depth estimation method that uses kernel regression
to refine output depth with normals. However, it is sensitive to noisy outputs
and is inherently incapable of handling cost volume inputs. Another line of
works [9, 15, 42] proposes the depth-normal consistency loss to enhance the
network’s perception of geometric cues. Unlike these methods, our method
leverages the normal to yield geometrically consistent costs in the 3D space,
yielding better utilization of adjacent costs.

Figure 2: Overview of our method. Given a reference image and a set of source
images, we use FPN to extract multi-scale features for cost volume
reconstruction. To conduct geometrically consistent aggregation within the
local window, we collect adjacent geometric cues and send them to the proposed
geometrically consistent propagation (GCP) module, which computes the
correspondence from the adjacent depth hypothesis space to the reference depth
space. The resulting costs are endowed with geometric consistency, which
facilitates better utilization of adjacent geometry and can be aggregated by
the convolution.

##  3 Methodology

Given a reference image $\mathbf{I}_{0}\in\mathbb{R}^{H\times W\times 3}$ and
$N$ source images $\\{\mathbf{I}_{i}\\}{{}_{i=1}^{N}}$, as well as camera
intrinsic $\\{\mathbf{K}_{i}\\}_{i=0}^{N}$ and extrinsic parameters
$\\{[\mathbf{R}_{0\rightarrow i};\mathbf{t}_{0\rightarrow i}]\\}_{i=1}^{N}$,
our goal is to estimate the depth map of $\mathbf{I}_{0}$ from multiple posed
images. Fig. 2 shows an overview of our method. We first utilize multi-scale
image features to build the cost volume (Sec. 3.1). We then introduce the
geometrically consistent aggregation scheme (Sec. 3.2), which consists of the
blocks in the depth network. We then investigate different choices for
obtaining surface normals (Sec. 3.3).

###  3.1 Cost Volume Construction

We first apply a Feature Pyramid Network [13] to extract multi-scale image
features
$\\{\mathbf{F}_{i}^{s}\\}_{i=0}^{N}\in\mathbb{R}^{\frac{H}{2^{s}}\times\frac{W}{2^{s}}\times
M}$, where $s$ is the scale factor. For simplicity, we omit the superscript of
$s$ below. To build the cost volume in each stage, we first sample depth
hypotheses $d$ for each pixel in a predefined depth range. Through
differentiable homography, we can compute the corresponding position
$\mathbf{p}^{\prime}$ of the reference image’s pixel $\mathbf{p}$ in the
source image,

$\mathbf{p}^{\prime}=\mathbf{K}_{i}[\mathbf{R}_{0\rightarrow
i}(\mathbf{K}_{0}^{-1}\mathbf{p}d)+\mathbf{t}_{0\rightarrow i}],$ (1)

where $\mathbf{R}$ and $\mathbf{t}$ denote the rotation and translation
parameters and $\mathbf{K}$ are the intrinsic matrix. Let
${\mathbf{F}(\mathbf{p})}$ represents the feature vector at pixel
$\mathbf{p}$, then the two-view feature correlation volume $\mathbf{V}$ at
pixel $\mathbf{p}$ can be represented as

$\mathbf{V}_{i}(\mathbf{p})=\mathbf{F}_{0}(\mathbf{p})\cdot\mathbf{F}_{i}({\mathbf{p}^{\prime}}),$
(2)

where $\cdot$ refers to the dot product. To aggregate multiple pair-wise cost
volumes, we utilize a shallow network [25] to learn the pixel-wise weight maps
$\mathbf{W}$. The weight computation takes place exclusively in the initial
stage, while weight maps for subsequent stages are derived through upsampling
from the previous stage. Then the multi-view aggregated cost volume
$\mathbf{C}$ can be represented as:

$\mathbf{C}=\frac{{\sum_{i=1}^{N}}{\mathbf{W}_{i}\odot\mathbf{V}_{i}}}{{\sum_{i=1}^{N}}{\mathbf{W}_{i}}}.$
(3)

###  3.2 Geometrically Consistent Aggregation

An essential idea of cost aggregation is to leverage neighboring information
to improve the discriminativeness of the cost volume, where the key is to find
the most relevant neighbors and effectively aggregate their matching costs. To
achieve this, typical convolution-based methods are limited to the size of the
convolution kernel (_e.g_., $3\times 3\times 3$), and geometric inconsistency
is very likely to happen in this local region due to non-constant depth
distributions within this kernel. It’s also computationally inefficient to
directly increase the kernel size to get improved performance.

In this paper, we observe in a small local region, many scenes can be
approximated with a plane, which frequently exists in real-world scenarios. To
this end, we propose to leverage this locally approximated planar structure to
guide the cost aggregation process in a geometrically consistent manner. There
exists an analytic relationship between the reference pixel’s depth and its
local neighbors, which can be leveraged to obtain more reliable cost
candidates. Specifically, for each reference pixel, we first collect the
geometric clues of its $k\times k$ spatial window to compute the
correspondences of the depth hypothesis. Depending on the corresponding
location, we propagate the adjacent costs to the reference pixel’s depth
space. Finally, we use a convolution layer to aggregate the propagated costs.

####  3.2.1 Local Geometric Clues Collection

We first collect local depth hypotheses and normal maps for each pixel within
a spatial window. Specifically, given the depth hypotheses of shape $L\times
H\times W$ and the normal map of shape $3\times H\times W$, where $L$ is the
depth hypothesis number and $H$, $W$ denotes the spatial dimension, we unfold
each pixel with a $k\times k$ spatial window, yielding local intermediate
depth hypotheses volume and normal map of shape $k^{2}\times L\times H\times
W$ and $k^{2}\times 3\times H\times W$, respectively. We then compute the
depth hypothesis correspondences based on these intermediate geometric clues.

####  3.2.2 Geometrically Consistent Propagation

To better aggregate the high-quality costs of the adjacent pixels, we align
the adjacent pixels’ depth hypothesis to the depth space of the reference
pixel. Based on depth correspondence, we perform geometrically consistent cost
propagation (GCP). Firstly, we introduce the depth relationship among pixels
within the same plane. Given a pixel’s image coordinates $(u,v)$ and depth
$d(u,v)$, its 3D point $X(u,v)$ in the camera coordinate system can be
represented as

$X(u,v)=\left[\begin{array}[]{ccc}x\\\ y\\\
z\end{array}\right]=\left[\begin{array}[]{ccc}\frac{u-c_{x}}{f_{x}}\\\
\frac{v-c_{y}}{f_{y}}\\\ {1}\end{array}\right]{d(u,v)},$ (4)

where $c_{x}$, $c_{y}$, $f_{x}$, and $f_{y}$ are the parameters of camera
intrinsic $\mathbf{K}$. For the given reference pixel $i$ and adjacent pixel
$j$, we model the relationship between $X(u_{i},v_{i})$ and $X(u_{j},v_{j})$
by leveraging local planar assumption and the surface normal $\mathbf{n}$.
They satisfy the equation of

${\mathbf{n}^{\top}}(X(u_{i},v_{i})-X(u_{j},v_{j}))=0.$ (5)

According to Eq. (4) and Eq. (5) , the depth relationship between the
reference pixel $i$ and the adjacent pixels $j$ can be represented as:

$\frac{d(u_{j},v_{j})}{d(u_{i},v_{i})}=\frac{{\mathbf{n}^{\top}}\left[\begin{array}[]{ccc}\frac{u_{i}-c_{x}}{f_{x}}&\frac{v_{i}-c_{y}}{f_{y}}&{1}\end{array}\right]^{\top}}{{\mathbf{n}^{\top}}\left[\begin{array}[]{ccc}\frac{u_{j}-c_{x}}{f_{x}}&\frac{v_{j}-c_{y}}{f_{y}}&{1}\end{array}\right]^{\top}}.$
(6)

We use $r_{ji}=\frac{d(u_{j},v_{j})}{d(u_{i},v_{i})}$ to denote the depth
ratio between $j$ and $i$, which describes the linear transformation of depth
within the plane. Based on this, we can compute the depth hypothesis
correspondences. Specifically, define $[d^{1}_{i},...,d^{L}_{i}]$ as the depth
hypothesis in the pixel $i$’s depth space, where $L$ refers to the number of
depth sampling levels. Each depth hypothesis is then mapped to pixel $j$’s
depth space through the depth ratio $r_{ji}$.

$[d^{1}_{i\rightarrow j},...,d^{L}_{i\rightarrow j}]=[r_{ji}\times
d^{1}_{i},...,r_{ji}\times d^{L}_{i}],$ (7)

where $d_{i\rightarrow j}$ represents the mapping depth of pixel $i$’s depth
hypothesis in pixel $j$’s depth space. We then propagate the matching cost of
pixel $j$ at the $d_{i\rightarrow j}$ to $d_{i}$. Let $\mathbf{C}_{j}$ denote
the cost for pixel $j$. The propagated matching cost $\mathbf{C}_{j\rightarrow
i}$ can be expressed as:

$\mathbf{C}_{j\rightarrow
i}(d^{0}_{i},...,d^{l}_{i})=\mathbf{C}_{j}(d^{0}_{i\rightarrow
j},...,d^{l}_{i\rightarrow j}).$ (8)

Since depth hypotheses are discretely sampled at regular depth intervals
within the depth range, we can conveniently use linear interpolation to
implement the above process. With the definition $d^{m}_{i\rightarrow
j}=d^{n}_{j}$, $\mathbf{C}_{j\rightarrow i}(d^{m}_{i})$ can be expressed as:

$\mathbf{C}_{j\rightarrow i}(d^{m}_{i})=(\mathbf{C}_{j}(d^{\lceil
n\rceil}_{j})-\mathbf{C}_{j}(d^{\lfloor n\rfloor}_{j}))\frac{n-\lfloor
n\rfloor}{\lceil n\rceil-\lfloor n\rfloor}.$ (9)

We refer to this process as geometrically consistent propagation from $j$ to
$i$. It can generate geometrically consistent cost candidates for each
reference pixel. Due to varying depth relationships between each pixel and its
adjacent pixels, cost propagation generates an intermediate cost of
$k^{2}M\times L\times H\times W$, where $M$ is the channel dimension.

####  3.2.3 Aggregating Propagated Costs

Since the intermediate costs include $k\times k$ spatial information in the
channel dimension, we thus aggregate the costs using convolutions with a
kernel size of $1\times 1\times k$ and an expanded channel dimension $k^{2}M$,
leading to the same parameters as the generic 3D convolutions with kernel size
$k\times k\times k$.

We encapsulate GCP and the convolution into one geometrically consistent
aggregation operator used to build the depth network. In particular, we still
keep the 3D U-Net architecture proposed by MVSNet [38], while replacing each
standard 3D convolution block with our proposed geometrically consistent
aggregation operator. For the upsampling layer in the U-Net structure, we use
the pixel shuffle to reorganize features and obtain a high-resolution cost
volume.

###  3.3 Extracting Normal Cues

Since our approach uses the surface normal for cost aggregation, in this
section, we study different methods for obtaining surface normals. We conduct
experiments to demonstrate the effectiveness of each method in Sec. 4.4.

Depth to normal. The surface normal can be directly computed from the
estimated depth. Since we use a three-stage cascade structure, we leverage the
depth map from the $g$ stage to generate the surface normal for the $g+1$
stage. The normal $\mathbf{n}$ can be computed [19] in closed form as:

$\mathbf{n}=\frac{(\mathbf{A}^{\mathsf{T}}\mathbf{A})^{-1}{\mathbf{A}^{\mathsf{T}}}{\mathbf{1}}}{\|{(\mathbf{A}^{\mathsf{T}}\mathbf{A})^{-1}{\mathbf{A}^{\mathsf{T}}}{\mathbf{1}}}\|},$
(10)

where $\mathbf{A}$ is a matrix composed of the coordinates of all pixels
within the local window. In addition to using estimated depth maps, we also
compute the GT normal from the GT depth maps following the same protocol and
use it to train our method for evaluating performances.

Cost to normal. In addition, inspired by [9], we use an additional network
branch to directly regress the normal map from the cost volume in each stage,
which is then used as a prior for geometrically consistent aggregation.

Off-the-shelf monocular surface normal. Monocular networks directly perceive
surface geometry from deep features and can estimate reasonable solutions in
regions with multi-view consistency ambiguities, which complements the task of
MVS. Therefore, we explore an existing monocular normal estimation network
Omnidata [5] to generate the surface normal. Since Omnidata is trained on low-
resolution images, its normal prediction might become unreliable when the
testing input resolution is increased. To tackle this, we adopt a divide-and-
conquer approach following MonoSDF [43] to generate high-resolution normal
cues. Specifically, we first divide the high-resolution image into multiple
overlapping patches. Surface normal estimation is then independently conducted
for each patch. Subsequently, the surface normal results are aligned and fused
to generate a high-resolution normal map.

###  3.4 Optimization

We treat the MVS task as a classification problem and employ the winner-takes-
all strategy to obtain the final depth map [39]. We use the cross-entropy loss
(Eq. 11) in each stage, which is applied to the probability volume $P$ and the
ground truth one-hot volume $P^{{}^{\prime}}$. Following [17], all depth out-
of-range will be masked during the training stage.

$\mathcal{L}=\sum_{i=1}^{d}{-}P^{{}^{\prime}}_{i}\log(P_{i}).$ (11)

(a) TransMVSNet [4] (b) GeoMVSNet [46] (c) Ours (d) GT

Figure 3: Comparison of reconstruction results. Our method reconstructs more
complete results in challenging areas.

##  4 Experiments

In this section, we evaluate our method on the DTU [1], ETH3D [22], and Tanks
and Temple [8] datasets, respectively. Furthermore, we conducted multiple
ablation experiments on the DTU dataset to validate the effectiveness of our
method.

###  4.1 Datasets

DTU [1] dataset comprises 128 scenes in controlled laboratory environments,
with models captured using structured light scanners. Each scene was scanned
from the same 49 or 64 camera positions under 7 different lighting conditions.
The official evaluation assesses the point cloud using distance metrics of
accuracy and completeness. BlendedMVS [40] is a large-scale MVS dataset that
consists of over 17,000 high-resolution images covering a variety of scenes,
including urban environments, architecture, sculptures, and small objects.
Tanks and Temples (TNT) [8] is a real-world dataset, divided into two sets,
including 8 scenes in the intermediate set and 6 scenes in the advanced set.
ETH3D [22] dataset consists of multiple indoor and outdoor scenes with large
viewpoint variations. The quality of point clouds on the ETH3D and TNT
datasets is measured using the percentage of precision and recall.

###  4.2 Implementation Details

##### Training

Following the data partitioning of MVSNet, we first train the model on the DTU
training set. Our network employs a three-stage cascade structure, with depth
sampling at 48, 32, and 8 in each stage and depth intervals of 4, 1, and 0.5,
respectively. We train our model with 5 input images, each having a resolution
of 512$\times$640\. The model is optimized using Adam for 12 epochs, starting
with an initial learning rate of 0.001 which is reduced by 0.5 after the 6 and
8 epochs. We then fine-tune the model on the BlendedMVS dataset with 9 images
at a resolution of 576$\times$768 for evaluation on Tanks and Temples and
ETH3D datasets. During fine-tuning, we reduce the depth sampling interval of
the last stage by half of its original value.

##### Evaluation

When testing on the DTU dataset, we use 5 images at a resolution of
864$\times$1152 as input and employ the depth map filtering method following
[46] to generate the final point cloud. For the tanks and temple dataset, we
carried out tests using 11 images with a resolution of 960$\times$1920\. In
terms of depth map fusion, we employ the widely adopted dynamic fusion
strategy [35]. Moreover, we conducted tests on the ETH3D dataset using images
with a size of 1152$\times$1536 and the depth map fusion strategy is
consistent with IterMVS [26].

(a) Family (b) Francis (c) Courtroom (d) Musume

Figure 4: Qualitative results on Tanks and Temples. Our method achieves
detailed and complete reconstructions across different scenes.

###  4.3 Benchmark Performance

##### Evaluation on DTU dataset.

Method | Acc. $\downarrow$ | Comp. $\downarrow$ | Overall$\downarrow$
---|---|---|---
Gipuma [6] | 0.283 | 0.873 | 0.578
COLMAP [21] | 0.400 | 0.664 | 0.532
NAPV-MVS [24] | 0.367 | 0.375 | 0.371
AA-RMVSNet [29] | 0.376 | 0.339 | 0.357
Vis-MVSNet [44] | 0.369 | 0.361 | 0.365
CasMVSNet [7] | 0.325 | 0.385 | 0.355
UniMVSNet [18] | 0.352 | 0.278 | 0.315
MVSTER [27] | 0.350 | 0.276 | 0.313
TransMVSNet [4] | 0.321 | 0.289 | 0.305
GbiNet* [17] | 0.314 | 0.295 | 0.305
RA-MVSNet [45] | 0.326 | 0.268 | 0.297
GeoMVSNet [46] | 0.331 | 0.259 | 0.295
ET-MVSNet [14] | 0.329 | 0.253 | 0.291
MVSformer [2] | 0.327 | 0.251 | 0.289
GoMVS | 0.347 | 0.227 | 0.287

Table 1: Quantitative results on DTU [1]. Our method achieves the best
completeness and overall score. Moreover, the completeness of our point cloud
outperforms previous methods by large margins.

We compare both traditional methods and deep learning-based approaches. The
quantitative evaluation results for point cloud reconstruction are shown in
Tab 1. Our method achieves SOTA completeness and overall performance. It is
worth noting that our method shows obvious improvement in completeness
compared to previous methods. This demonstrates that our method can better use
adjacent costs to propagate local geometries, resulting in a more complete
reconstruction. Fig. 3 shows a comparison of our point cloud results with
previous SOTA methods. We have more detailed and complete reconstructions in
the challenge areas.

Methods | Intermediate | Advanced
---|---|---
Mean$\uparrow$ | Fam. | Fra. | Hor. | Lig. | M60 | Pan. | Pla. | Tra. | Mean$\uparrow$ | Aud. | Bal. | Cou. | Mus. | Pal. | Tem.
COLMAP [21] | 42.14 | 50.41 | 22.25 | 26.63 | 56.43 | 44.83 | 46.97 | 48.53 | 42.04 | 27.24 | 16.02 | 25.23 | 34.70 | 41.51 | 18.05 | 27.94
CasMVSNet [7] | 56.84 | 76.37 | 58.45 | 46.26 | 55.81 | 56.11 | 54.06 | 58.18 | 49.51 | 31.12 | 19.81 | 38.46 | 29.10 | 43.87 | 27.36 | 28.11
Vis-MVSNet [44] | 60.03 | 77.40 | 60.23 | 47.07 | 63.44 | 62.21 | 57.28 | 60.54 | 52.07 | 33.78 | 20.79 | 38.77 | 32.45 | 44.20 | 28.73 | 37.70
GBiNet [17] | 61.42 | 79.77 | 67.69 | 51.81 | 61.25 | 60.37 | 55.87 | 60.67 | 53.89 | 37.32 | 29.77 | 42.12 | 36.30 | 47.69 | 31.11 | 36.93
EPP-MVSNet [16] | 61.68 | 77.86 | 60.54 | 52.96 | 62.33 | 61.69 | 60.34 | 62.44 | 55.30 | 35.72 | 21.28 | 39.74 | 35.34 | 49.21 | 30.00 | 38.75
TransMVSNet [4] | 63.52 | 80.92 | 65.83 | 56.94 | 62.54 | 63.06 | 60.00 | 60.20 | 58.67 | 37.00 | 24.84 | 44.59 | 34.77 | 46.49 | 34.69 | 36.62
UniMVSNet [18] | 64.36 | 81.20 | 66.43 | 53.11 | 64.36 | 66.09 | 64.84 | 62.23 | 57.53 | 38.96 | 28.33 | 44.36 | 39.74 | 52.89 | 33.80 | 34.63
D-MVSNet [41] | 64.66 | 81.27 | 67.54 | 59.10 | 63.12 | 64.64 | 64.80 | 59.83 | 56.97 | 41.17 | 30.08 | 46.10 | 40.65 | 53.53 | 35.08 | 41.60
ET-MVSNet [14] | 65.49 | 81.65 | 68.79 | 59.46 | 65.72 | 64.22 | 64.03 | 61.23 | 58.79 | 40.41 | 28.86 | 45.18 | 38.66 | 51.10 | 35.39 | 43.23
RA-MVSNet [45] | 65.72 | 82.44 | 66.61 | 58.40 | 64.78 | 67.14 | 65.60 | 62.74 | 58.08 | 39.93 | 29.17 | 46.05 | 40.23 | 53.22 | 34.62 | 36.30
GeoMVSNet [46] | 65.89 | 81.64 | 67.53 | 55.78 | 68.02 | 65.49 | 67.19 | 63.27 | 58.22 | 41.52 | 30.23 | 46.53 | 39.98 |  53.05 | 35.98 | 43.34
MVSFormer [2] | 66.37 | 82.06 | 69.34 | 60.49 | 68.61 | 65.67 | 64.08 | 61.23 | 59.33 | 40.87 | 28.22 | 46.75 | 39.30 |  52.88 | 35.16 | 42.95
GoMVS | 66.44 | 82.68 | 69.23 | 69.19 | 63.56 | 65.13 | 62.10 | 58.81 | 60.80 | 43.07 | 35.52 | 47.15 | 42.52 | 52.08 | 36.34 | 44.82

Table 2: Quantitative results of F-score on Tanks and Temples benchmark. Our
method achieves the best F-score on both the “Intermediate” and the
challenging “Advanced” set. Note that our method ranks 1st on the official TNT
Advanced Benchmark.

##### Evaluation on Tanks and Temples dataset.

We validated the generalization of our model on the Tanks and Temples dataset,
and the quantitative results are shown in Table . We achieved the best
performance on both the intermediate and advanced sets. Moreover, we ranked
$1st$ among all submitted results on the advanced set of the TNT benchmark,
which contains more complex scenes. It demonstrates the strong robustness and
generalization ability of our method. Fig. 4 shows point cloud results on
intermediate and advanced sets. Our method achieves detailed and complete
reconstructions across different indoor and outdoor scenes.

Method | Acc. $\downarrow$ | Comp. $\downarrow$ | Overall $\downarrow$
---|---|---|---
Standard 3D-convolution [38] | 0.365 | 0.265 | 0.315
Sparse convolution [37] | 0.354 | 0.268 | 0.311
Spatial deformable aggregation [25] | 0.363 | 0.257 | 0.310
Depth kernel regression [19] | 0.369 | 0.262 | 0.316
Ours (GCA) | 0.347 | 0.227 | 0.287

Table 3: Comparison with different aggregation methods. Our method
significantly outperforms previous cost volume aggregation methods.

##### Evaluation on ETH3D dataset.

The ETH3D dataset contains many challenging scenes, including scenes with
textureless areas and large viewpoint variations. We compare our methods with
previous methods and results are shown in Tab. 4. Our method achieves the best
performance on both the validation set and the test split. In particular, it
outperforms previous SOTA by a significant margin on the test split,
demonstrating its generalization ability over existing methods.

Methods |  Training |  Test
---|---|---
Acc.$\uparrow$ |  Comp. $\uparrow$ |  F-score$\uparrow$ |  Acc.$\uparrow$ |  Comp. $\uparrow$ |  F-score$\uparrow$
COLMAP [21] |  91.85 |  55.13 |  67.66 |  91.97 |  62.98 |  73.01
ACMM [32] |  90.67 |  70.42 |  78.86 |  90.65 |  74.34 |  80.78
IterMVS [26] |  73.62 |  61.87 |  66.36 |  76.91 |  72.65 |  74.29
GBi-Net [17] |  73.17 |  69.21 |  70.78 |  80.02 |  75.65 |  78.40
MVSTER [27] |  76.92 |  68.08 |  72.06 |  77.09 |  82.47 |  79.01
PVSNet [34] |  83.00 |  71.76 |  76.57 |  81.55 |  83.97 |  82.62
EPP-MVSNet [16] |  82.76 |  67.58 |  74.00 |  85.47 |  81.79 |  83.40
Vis-MVSNet [44] |  83.32 |  65.53 |  72.77 |  86.86 |  80.92 |  83.46
EPNet [23] |  79.36 |  79.28 |  79.08 |  80.37 |  87.84 |  83.72
GoMVS |  81.22 |  77.65 |  79.16 |  85.50 |  86.85 |  85.91

Table 4: Quantitative results on ETH3D dataset. We show comparisons of
reconstructed point clouds using percentage metric (%) at a threshold of 2cm.
Our approach achieves the best performance with notable margins.

###  4.4 Ablation Study

##### Comparison with different aggregation methods.

To verify the effectiveness of utilizing adjacent geometry, we compare
different cost aggregation and depth aggregation methods, and the results are
shown in Tab. 3. Regarding the cost aggregation methods, Sparse convolution
[37] aggregates the cost at the same depth without fully considering the depth
geometry, resulting in certain improvements in performance compared with the
baseline. PatchMatchNet [25] utilizes deformable convolutions to gather
spatial matching costs and aggregate them using a lightweight 3D CNN. We
replace the aggregation network with a 3D U-Net to ensure a fair comparison
with the same parameter scale. PatchmatchNet heavily relies on network
capabilities and does not guarantee geometric plausibility from the selected
costs. As a result, it brings limited performance improvements (row #3).

Additionally, for the depth aggregation method, we refine the depth map on the
baseline method by incorporating depth kernel regression proposed by GeoNet
[19]. Using normal similarity to compute depth aggregation weights is prone to
the influence of normal noise and cannot effectively utilize the abundant
geometric information in the cost volume. This leads to a decline in the
accuracy of the final point cloud (row #4). We utilize normal priors to guide
cost aggregation, alleviating the challenge of geometric inconsistency and
achieving the best performance among all aggregation methods.

##### Comparison with different depth receptive fields.

Intuitively, 3D convolutions with larger receptive fields in the depth
dimension can alleviate the cost inconsistency in the local range, by
resorting to wider areas. Therefore, we compare our approach with variants
directly expanding the depth receptive field. We keep the $3\times 3$ spatial
window size at each 3D convolution layer and experiment with kernel sizes of
3, 5, and 7 in the depth dimension on the baseline method. The quantitative
results are shown in Tab. 5, we find that increasing the receptive field in
the depth dimension leads to some certain improvement. However, due to the
lack of geometric awareness, its performance is saturated when the dimension
expands to a certain kernel size. In contrast, we use surface normal to
geometrically guide the cost aggregation process. With a kernel size of only
3, our method achieves the best performance, outperforming other alternatives
by clear margins.

##### Evaluation of different normal cues.

Since the surface normal is important for guiding geometrically consistent
aggregation, we further evaluate the effectiveness of different normal cues in
Tab. 6. We first train and evaluate our method using the GT normal, which sets
an upper bound for our method. As shown in the last row, it significantly
improves the performance of point clouds, validating our method’s
effectiveness when using high-quality normal inputs. We further train and
evaluate our method using depth-computed normals [19] or cost-computed normals
[9], the results are suboptimal as they essentially rely on the quality of
input depth, which can degrade in challenging areas. Though lacking multi-view
consistency, monocular normals do not collapse in challenging geometric
estimation regions of the cost volume. This reveals a nice property for
monocular estimations. In addition to the DTU dataset, we also observe notable
improvement using monocular surface normals on other benchmarks.

Aggregation kernel $(d\times h\times w)$ | Acc. $\downarrow$ | Comp. $\downarrow$ | Overall $\downarrow$
---|---|---|---
Standard Conv3D$(3\times 3\times 3)$ | 0.365 | 0.265 | 0.315
Standard Conv3D $(5\times 3\times 3)$ | 0.352 | 0.260 | 0.306
Standard Conv3D $(7\times 3\times 3)$ | 0.352 | 0.258 | 0.305
Proposed GCA $(3\times 3\times 3)$ | 0.347 | 0.227 | 0.287

Table 5: Evaluation of aggregation receptive fields. Directly expanding
receptive fields along the depth dimension yields limited improvement and is
easily saturated. In contrast, our method achieves the best performance with a
kernel size of 3.

Method | Acc. $\downarrow$ | Comp. $\downarrow$ | Overall $\downarrow$
---|---|---|---
Ours + Depth-to-normal [19] | 0.352 | 0.242 | 0.297
Ours + Cost-to-normal [9] | 0.358 | 0.241 | 0.300
Ours + Mono-normal [5] | 0.347 | 0.227 | 0.287
Ours + GT normal | 0.275 | 0.221 | 0.248

Table 6: Evaluation of different normal cues. Our method with GT normal
demonstrates remarkable performance (0.248). Among all estimated normals, the
off-the-shelf monocular normal has the best performance.

##  5 Conclusion

In this paper, we propose GoMVS, which aggregates locally consistent
geometries to better utilize adjacent geometry. By leveraging local smoothness
in conjunction with surface normal, we propose geometrically consistent
aggregation. it computes the correspondence from the adjacent depth hypotheses
space to the reference depth space and propagates cost accordingly.
Furthermore, we investigate different choices for generating normal priors and
find that monocular cues effectively complement the MVS network. Our method
achieves state-of-the-art performance on the DTU, Tanks and Temples, and ETH3D
datasets.

Acknowledgements Y. Zhang was supported by NSFC (No.U19B2037) and the Natural
Science Basic Research Program of Shaanxi (No.2021JCW-03). Y. Zhu was
supported by NSFC (No.61901384).

## References

  * Aanæs et al. [2016] Henrik Aanæs, Rasmus Ramsbøl Jensen, George Vogiatzis, Engin Tola, and Anders Bjorholm Dahl.  Large-scale data for multiple-view stereopsis.  _International Journal of Computer Vision_ , 120:153–168, 2016.
  * Cao et al. [2022] Chenjie Cao, Xinlin Ren, and Yanwei Fu.  Mvsformer: Learning robust image representations via transformers and temperature-based depth for multi-view stereo.  _arXiv preprint arXiv:2208.02541_ , 2022.
  * Cheng et al. [2020] Shuo Cheng, Zexiang Xu, Shilin Zhu, Zhuwen Li, Li Erran Li, Ravi Ramamoorthi, and Hao Su.  Deep stereo using adaptive thin volume representation with uncertainty awareness.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 2524–2534, 2020.
  * Ding et al. [2022] Yikang Ding, Wentao Yuan, Qingtian Zhu, Haotian Zhang, Xiangyue Liu, Yuanjiang Wang, and Xiao Liu.  Transmvsnet: Global context-aware multi-view stereo network with transformers.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 8585–8594, 2022.
  * Eftekhar et al. [2021] Ainaz Eftekhar, Alexander Sax, Jitendra Malik, and Amir Zamir.  Omnidata: A scalable pipeline for making multi-task mid-level vision datasets from 3d scans.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 10786–10796, 2021.
  * Galliani et al. [2015] Silvano Galliani, Katrin Lasinger, and Konrad Schindler.  Massively parallel multiview stereopsis by surface normal diffusion.  In _Proceedings of the IEEE International Conference on Computer Vision_ , pages 873–881, 2015.
  * Gu et al. [2020] Xiaodong Gu, Zhiwen Fan, Siyu Zhu, Zuozhuo Dai, Feitong Tan, and Ping Tan.  Cascade cost volume for high-resolution multi-view stereo and stereo matching.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 2495–2504, 2020.
  * Knapitsch et al. [2017] Arno Knapitsch, Jaesik Park, Qian-Yi Zhou, and Vladlen Koltun.  Tanks and temples: Benchmarking large-scale scene reconstruction.  _ACM Transactions on Graphics (ToG)_ , 36(4):1–13, 2017.
  * Kusupati et al. [2020] Uday Kusupati, Shuo Cheng, Rui Chen, and Hao Su.  Normal assisted stereo depth estimation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 2189–2199, 2020.
  * Li et al. [2023a] Jingliang Li, Zhengda Lu, Yiqun Wang, Jun Xiao, and Ying Wang.  Nr-mvsnet: Learning multi-view stereo based on normal consistency and depth refinement.  _IEEE Transactions on Image Processing_ , 2023a.
  * Li et al. [2023b] Rui Li, Dong Gong, Wei Yin, Hao Chen, Yu Zhu, Kaixuan Wang, Xiaozhi Chen, Jinqiu Sun, and Yanning Zhang.  Learning to fuse monocular and multi-view cues for multi-frame depth estimation in dynamic scenes.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 21539–21548, 2023b.
  * Liao et al. [2022] Jinli Liao, Yikang Ding, Yoli Shavit, Dihe Huang, Shihao Ren, Jia Guo, Wensen Feng, and Kai Zhang.  Wt-mvsnet: window-based transformers for multi-view stereo.  _Advances in Neural Information Processing Systems_ , 35:8564–8576, 2022.
  * Lin et al. [2017] Tsung-Yi Lin, Piotr Dollár, Ross Girshick, Kaiming He, Bharath Hariharan, and Serge Belongie.  Feature pyramid networks for object detection.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 2117–2125, 2017.
  * Liu et al. [2023] Tianqi Liu, Xinyi Ye, Weiyue Zhao, Zhiyu Pan, Min Shi, and Zhiguo Cao.  When epipolar constraint meets non-local operators in multi-view stereo.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 18088–18097, 2023.
  * Long et al. [2020] Xiaoxiao Long, Lingjie Liu, Christian Theobalt, and Wenping Wang.  Occlusion-aware depth estimation with adaptive normal constraints.  In _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part IX 16_ , pages 640–657. Springer, 2020.
  * Ma et al. [2021] Xinjun Ma, Yue Gong, Qirui Wang, Jingwei Huang, Lei Chen, and Fan Yu.  Epp-mvsnet: Epipolar-assembling based depth prediction for multi-view stereo.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 5732–5740, 2021.
  * Mi et al. [2022] Zhenxing Mi, Chang Di, and Dan Xu.  Generalized binary search network for highly-efficient multi-view stereo.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 12991–13000, 2022.
  * Peng et al. [2022] Rui Peng, Rongjie Wang, Zhenyu Wang, Yawen Lai, and Ronggang Wang.  Rethinking depth estimation for multi-view stereo: A unified representation.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 8645–8654, 2022.
  * Qi et al. [2018] Xiaojuan Qi, Renjie Liao, Zhengzhe Liu, Raquel Urtasun, and Jiaya Jia.  Geonet: Geometric neural network for joint depth and surface normal estimation.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 283–291, 2018.
  * Ren et al. [2023] Chunlin Ren, Qingshan Xu, Shikun Zhang, and Jiaqi Yang.  Hierarchical prior mining for non-local multi-view stereo.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 3611–3620, 2023.
  * Schönberger et al. [2016] Johannes L Schönberger, Enliang Zheng, Jan-Michael Frahm, and Marc Pollefeys.  Pixelwise view selection for unstructured multi-view stereo.  In _Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11-14, 2016, Proceedings, Part III 14_ , pages 501–518. Springer, 2016.
  * Schöps et al. [2017] Thomas Schöps, Johannes L. Schönberger, Silvano Galliani, Torsten Sattler, Konrad Schindler, Marc Pollefeys, and Andreas Geiger.  A multi-view stereo benchmark with high-resolution images and multi-camera videos.  In _Conference on Computer Vision and Pattern Recognition (CVPR)_ , 2017.
  * Su and Tao [2023] Wanjuan Su and Wenbing Tao.  Efficient edge-preserving multi-view stereo network for depth estimation.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , pages 2348–2356, 2023.
  * Tong et al. [2022] Wei Tong, Xiaorong Guan, Jian Kang, Poly ZH Sun, Rob Law, Pedram Ghamisi, and Edmond Q Wu.  Normal assisted pixel-visibility learning with cost aggregation for multiview stereo.  _IEEE Transactions on Intelligent Transportation Systems_ , 23(12):24686–24697, 2022.
  * Wang et al. [2021] Fangjinhua Wang, Silvano Galliani, Christoph Vogel, Pablo Speciale, and Marc Pollefeys.  Patchmatchnet: Learned multi-view patchmatch stereo.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 14194–14203, 2021.
  * Wang et al. [2022a] Fangjinhua Wang, Silvano Galliani, Christoph Vogel, and Marc Pollefeys.  Itermvs: Iterative probability estimation for efficient multi-view stereo.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 8606–8615, 2022a.
  * Wang et al. [2022b] Xiaofeng Wang, Zheng Zhu, Guan Huang, Fangbo Qin, Yun Ye, Yijia He, Xu Chi, and Xingang Wang.  Mvster: epipolar transformer for efficient multi-view stereo.  In _Computer Vision–ECCV 2022: 17th European Conference, Tel Aviv, Israel, October 23–27, 2022, Proceedings, Part XXXI_ , pages 573–591. Springer, 2022b.
  * Wang et al. [2022c] Yun Wang, Longguang Wang, Hanyun Wang, and Yulan Guo.  Spnet: Learning stereo matching with slanted plane aggregation.  _IEEE Robotics and Automation Letters_ , 7(3):6258–6265, 2022c.
  * Wei et al. [2021] Zizhuang Wei, Qingtian Zhu, Chen Min, Yisong Chen, and Guoping Wang.  Aa-rmvsnet: Adaptive aggregation recurrent multi-view stereo network.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ , pages 6187–6196, 2021.
  * Wu et al. [2024] Jiang Wu, Rui Li, Yu Zhu, Wenxun Zhao, Jinqiu Sun, and Yanning Zhang.  Boosting multi-view stereo with late cost aggregation.  _arXiv preprint arXiv:2401.11751_ , 2024.
  * Xu and Zhang [2020] Haofei Xu and Juyong Zhang.  Aanet: Adaptive aggregation network for efficient stereo matching.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 1959–1968, 2020.
  * Xu and Tao [2019] Qingshan Xu and Wenbing Tao.  Multi-scale geometric consistency guided multi-view stereo.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 5483–5492, 2019.
  * Xu and Tao [2020] Qingshan Xu and Wenbing Tao.  Planar prior assisted patchmatch multi-view stereo.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , pages 12516–12523, 2020.
  * Xu et al. [2022] Qingshan Xu, Wanjuan Su, Yuhang Qi, Wenbing Tao, and Marc Pollefeys.  Learning inverse depth regression for pixelwise visibility-aware multi-view stereo networks.  _International Journal of Computer Vision_ , 130(8):2040–2059, 2022.
  * Yan et al. [2020] Jianfeng Yan, Zizhuang Wei, Hongwei Yi, Mingyu Ding, Runze Zhang, Yisong Chen, Guoping Wang, and Yu-Wing Tai.  Dense hybrid recurrent multi-view stereo net with dynamic consistency checking.  In _Computer Vision–ECCV 2020: 16th European Conference, Glasgow, UK, August 23–28, 2020, Proceedings, Part IV_ , pages 674–689. Springer, 2020.
  * Yang et al. [2020] Jiayu Yang, Wei Mao, Jose M Alvarez, and Miaomiao Liu.  Cost volume pyramid based depth inference for multi-view stereo.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 4877–4886, 2020.
  * Yang et al. [2022] Jiayu Yang, Jose M Alvarez, and Miaomiao Liu.  Non-parametric depth distribution modelling based depth inference for multi-view stereo.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 8626–8634, 2022.
  * Yao et al. [2018] Yao Yao, Zixin Luo, Shiwei Li, Tian Fang, and Long Quan.  Mvsnet: Depth inference for unstructured multi-view stereo.  In _Proceedings of the European conference on computer vision (ECCV)_ , pages 767–783, 2018.
  * Yao et al. [2019] Yao Yao, Zixin Luo, Shiwei Li, Tianwei Shen, Tian Fang, and Long Quan.  Recurrent mvsnet for high-resolution multi-view stereo depth inference.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pages 5525–5534, 2019.
  * Yao et al. [2020] Yao Yao, Zixin Luo, Shiwei Li, Jingyang Zhang, Yufan Ren, Lei Zhou, Tian Fang, and Long Quan.  Blendedmvs: A large-scale dataset for generalized multi-view stereo networks.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 1790–1799, 2020.
  * Ye et al. [2023] Xinyi Ye, Weiyue Zhao, Tianqi Liu, Zihao Huang, Zhiguo Cao, and Xin Li.  Constraining depth map geometry for multi-view stereo: A dual-depth approach with saddle-shaped depth cells.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision (ICCV)_ , pages 17661–17670, 2023.
  * Yin et al. [2019] Wei Yin, Yifan Liu, Chunhua Shen, and Youliang Yan.  Enforcing geometric constraints of virtual normal for depth prediction.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pages 5684–5693, 2019.
  * Yu et al. [2022] Zehao Yu, Songyou Peng, Michael Niemeyer, Torsten Sattler, and Andreas Geiger.  Monosdf: Exploring monocular geometric cues for neural implicit surface reconstruction.  _Advances in neural information processing systems_ , 35:25018–25032, 2022.
  * Zhang et al. [2023a] Jingyang Zhang, Shiwei Li, Zixin Luo, Tian Fang, and Yao Yao.  Vis-mvsnet: Visibility-aware multi-view stereo network.  _International Journal of Computer Vision_ , 131(1):199–214, 2023a.
  * Zhang et al. [2023b] Yisu Zhang, Jianke Zhu, and Lixiang Lin.  Multi-view stereo representation revist: Region-aware mvsnet.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 17376–17385, 2023b.
  * Zhang et al. [2023c] Zhe Zhang, Rui Peng, Yuxi Hu, and Ronggang Wang.  Geomvsnet: Learning multi-view stereo with geometry perception.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 21508–21518, 2023c.


# Multiple Kernel Learning: A Unifying Probabilistic Viewpoint

Hannes Nickisch<EMAIL_ADDRESS>
Max Planck Institute for Intelligent Systems, Spemannstraße 38, 72076
Tübingen, Germany
Matthias Seeger<EMAIL_ADDRESS>
Ecole Polytechnique Fédérale de Lausanne, INJ 339, Station 14, 1015 Lausanne,
Switzerland

###### Abstract

We present a probabilistic viewpoint to multiple kernel learning unifying
well-known regularised risk approaches and recent advances in approximate
Bayesian inference relaxations. The framework proposes a general objective
function suitable for regression, robust regression and classification that is
lower bound of the marginal likelihood and contains many regularised risk
approaches as special cases. Furthermore, we derive an efficient and provably
convergent optimisation algorithm.

Keywords: Multiple kernel learning, approximate Bayesian inference, double
loop algorithms, Gaussian processes

##  1 Introduction

Nonparametric kernel methods, cornerstones of machine learning today, can be
seen from different angles: as regularised risk minimisation in function
spaces (Schölkopf and Smola, 2002), or as probabilistic Gaussian process
methods (Rasmussen and Williams, 2006). In these techniques, the kernel (or
equivalently covariance) function encodes interpolation characteristics from
observed to unseen points, and two basic statistical problems have to be
mastered. First, a latent function must be predicted which fits data well, yet
is as smooth as possible given the fixed kernel. Second, the kernel function
parameters have to be learned as well, to best support predictions which are
of primary interest. While the first problem is simpler and has been addressed
much more frequently so far, the central role of learning the covariance
function is well acknowledged, and a substantial number of methods for
“learning the kernel”, “multiple kernel learning”, or “evidence maximisation”
are available now. However, much of this work has firmly been associated with
one of the “camps” (referred to as _regularised risk_ and _probabilistic_ in
the sequel) with surprisingly little crosstalk or acknowledgments of prior
work across this boundary. In this paper, we clarify the relationship between
major regularised risk and probabilistic kernel learning techniques precisely,
pointing out advantages and pitfalls of either, as well as algorithmic
similarities leading to novel powerful algorithms.

We develop a common analytical and algorithmical framework encompassing
approaches from both camps and provide clear insights into the optimisation
structure. Even though, most of the optimisation is non convex, we show how to
operate a provably convergent “almost Newton” method nevertheless. Each step
is not much more expensive than a gradient based approach. Also, we do not
require any foreign optimisation code to be available. Our framework unifies
kernel learning for regression, robust regression and classification.

The paper is structured as follows: In section 2, we introduce the regularised
risk and the probabilistic view of kernel learning. In increasing order of
generality, we explain multiple kernel learning (MKL, section 2.1), maximum a
posteriori estimation (MAP, section 2.2) and marginal likelihood maximisation
(MLM, section 2.3). A taxonomy of the mutual relations between the approaches
and important special cases is given in section 2.4. Section 3 introduces a
general optimisation scheme and section 4 draws a conclusion.

##  2 Kernel Methods and Kernel Learning

Kernel-based algorithms come in many shapes, however, the primary goal is –
based on training data $\\{(\mathbf{x}_{i},y_{i})\,|\,i=1..n\\}$,
$\mathbf{x}_{i}\in{\cal X}$, $y_{i}\in\mathcal{Y}$ and a parametrised kernel
function $k_{\bm{\theta}}(\mathbf{x},\mathbf{x}^{\prime})$ – to predict the
output $y_{*}$ for unseen inputs $\mathbf{x}_{*}$. Often, linear
parametrisations
$k_{\bm{\theta}}(\mathbf{x},\mathbf{x}^{\prime})=\sum_{m=1}^{M}\theta_{m}k_{m}(\mathbf{x},\mathbf{x}^{\prime})$
are used, where the $k_{m}$ are fixed positive definite functions, and
$\bm{\theta}\succeq\mathbf{0}$. Learning the kernel means finding
$\bm{\theta}$ to best support this goal. In general, kernel methods employ a
postulated latent function $u:{\cal X}\to\mathbb{R}$ whose smoothness is
controlled via the function space squared norm
$\|u(\cdot)\|_{k_{\bm{\theta}}}^{2}$. Most often, smoothness is traded against
data fit, either enforced by a _loss function_ $\ell(y_{i},u(\mathbf{x}_{i}))$
or modeled by a _likelihood_ $\mathbb{P}(y_{i}|u_{i})$. Let us define kernel
matrices
$\mathbf{K}_{\bm{\theta}}:=[k_{\bm{\theta}}(\mathbf{x}_{i},\mathbf{x}_{j})]_{ij}$,
and $\mathbf{K}_{m}:=[k_{m}(\mathbf{x}_{i},\mathbf{x}_{j})]_{ij}$ in
$\mathbb{R}^{n\times n}$ and the vectors
$\mathbf{y}:=[y_{i}]_{i}\in\mathcal{Y}^{n}$,
$\mathbf{u}:=[u(\mathbf{x}_{i})]_{i}\in\mathbb{R}^{n}$ collecting outputs and
latent function values, respectively.

The _regularised risk_ route to kernel prediction, which is followed by any
support vector machine (SVM) or ridge regression technique, yields
$\|u(\cdot)\|_{k_{\bm{\theta}}}^{2}+\frac{C}{n}\sum_{i=1}^{n}\ell(y_{i},u_{i})$
as criterion, enforcing smoothness of $u(\cdot)$ as well as good data fit
through th _e_ loss function $\frac{C}{n}\ell(y_{i},u(\mathbf{x}_{i}))$. By
the representer theorem, the minimiser can be written as
$u(\cdot)=\sum_{i}\alpha_{i}k_{\bm{\theta}}(\cdot,\mathbf{x}_{i})$, so that
$\|u(\cdot)\|_{k_{\bm{\theta}}}^{2}=\bm{\alpha}^{\top}\mathbf{K}_{\bm{\theta}}\bm{\alpha}$
(Schölkopf and Smola, 2002). As
$\mathbf{u}=\mathbf{K}_{\bm{\theta}}\bm{\alpha}$, the regularised risk problem
is given by

$\min_{\mathbf{u}}\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}+\frac{C}{n}\sum_{i=1}^{n}\ell(y_{i},u_{i}).$
(1)

A _probabilistic_ viewpoint of the same setting is based on the notion of a
Gaussian process (GP) (Rasmussen and Williams, 2006): a Gaussian random
function $u(\cdot)$ with mean function
$\mathbb{E}[u(\mathbf{x})]=m(\mathbf{x})\equiv 0$ and covariance function
$\mathbb{V}[u(\mathbf{x}),u(\mathbf{x}^{\prime})]=\mathbb{E}[u(\mathbf{x})u(\mathbf{x}^{\prime})]=k_{\bm{\theta}}(\mathbf{x},\mathbf{x}^{\prime})$.
In practice, we only use finite-dimensional snapshots of the process
$u(\cdot)$: for example,
$\mathbb{P}(\mathbf{u};\bm{\theta})=\mathcal{N}(\mathbf{u}|\mathbf{0},\mathbf{K}_{\bm{\theta}})$,
a zero-mean joint Gaussian with covariance matrix $\mathbf{K}_{\bm{\theta}}$.
We adopt this GP as prior distribution over $u(\cdot)$, estimating the latent
function as maximum of the posterior process
$\mathbb{P}(u(\cdot)|\mathbf{y};\bm{\theta})\propto\mathbb{P}(\mathbf{y}|\mathbf{u})\mathbb{P}(u(\cdot);\bm{\theta})$.
Since the likelihood depends on $u(\cdot)$ only through the finite subset
$\\{u(\mathbf{x}_{i})\\}$, the posterior process has a finite-dimensional
representation:
$\mathbb{P}(u(\cdot)|\mathbf{y},\mathbf{u})=\mathbb{P}(u(\cdot)|\mathbf{u})$,
so that
$\mathbb{P}(u(\cdot)|\mathbf{y};\bm{\theta})=\int\mathbb{P}(u(\cdot)|\mathbf{u})\mathbb{P}(\mathbf{u}|\mathbf{y};\bm{\theta})\text{d}\mathbf{u}$
is specified by the joint distribution
$\mathbb{P}(\mathbf{u}|\mathbf{y};\bm{\theta})$, a probabilistic equivalent of
the representer theorem. Kernel prediction amounts to _maximum a posteriori_
(MAP) estimation

$\max\nolimits_{\mathbf{u}}\mathbb{P}(\mathbf{u}|\mathbf{y};\bm{\theta})\equiv\max\nolimits_{\mathbf{u}}\mathbb{P}(\mathbf{u};\bm{\theta})\mathbb{P}(\mathbf{y}|\mathbf{u})\equiv\min\nolimits_{\mathbf{u}}\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}-2\ln\mathbb{P}(\mathbf{y}|\mathbf{u})+\ln|\mathbf{K}_{\bm{\theta}}|,$
(2)

ignoring an additive constant. Minimising equations (1) and (2) for any fixed
kernel matrix $\mathbf{K}$ gives the same minimiser $\hat{\mathbf{u}}$ and
prediction
$u(\mathbf{x}_{*})=\hat{\mathbf{u}}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}[k_{\bm{\theta}}(\mathbf{x}_{i},\mathbf{x}_{*})]_{i}$.

The correspondence between likelihood and loss function bridges probabilistic
and regularised risk techniques. More specifically, any likelihood
$\mathbb{P}(\mathbf{y}|\mathbf{u})$ induces a loss function
$\ell(\mathbf{y},\mathbf{u})$ via

$-2\ln\mathbb{P}(\mathbf{y}|\mathbf{u})=-2\sum_{i}\ln\mathbb{P}(y_{i}|u_{i})\rightsquigarrow\frac{C}{n}\sum_{i=1}^{n}\ell(y_{i},u_{i})=\ell(\mathbf{y},\mathbf{u}),$

however some loss functions cannot be interpreted as a negative log likelihood
as shown in table (2) and as discussed for the SVM by Sollich (2000). If, the
likelihood is a _log-concave_ function of $\mathbf{u}$, it corresponds to a
convex loss function (Boyd and Vandenberghe, 2002, Sect. 3.5.1). Common loss
functions and likelihoods for classification $\mathcal{Y}=\\{\pm 1\\}$ and
regression $\mathcal{Y}=\mathbb{R}$ are listed in table (2).

$\mathcal{Y}$ | Loss function | $\ell(y_{i},u_{i})$ | $\mathbb{P}(y_{i}|u_{i})$ | Likelihood
---|---|---|---|---
$\\{\pm 1\\}$ | SVM Hinge loss | $\max(0,1-y_{i}u_{i})$ | $\nexists$
$\\{\pm 1\\}$ | Log loss | $\ln(\exp(-y_{i}u_{i})+1)$ | $1/(\exp(-\tau y_{i}u_{i})+1)$ | Logistic
$\mathbb{R}$ | SVM $\epsilon$-insensitive loss | $\max(0,|y_{i}-u_{i}|/\epsilon-1)$ | $\nexists$
$\mathbb{R}$ | Quadratic loss | $(y_{i}-u_{i})^{2}$ | $\mathcal{N}(y_{i}|u_{i},\sigma^{2})$ | Gaussian
$\mathbb{R}$ | Linear loss | $|y_{i}-u_{i}|$ | $\mathcal{L}(y_{i}|u_{i},\tau)$ | Laplace

Table 1: Relations between loss functions and likelihoods

In the following, we discuss several approaches to learn the kernel parameters
$\bm{\theta}$ and show how all of them can be understood as instances of or
approximations to Bayesian evidence maximisation. Although the exposition MKL
section 2.1 and MAP section 2.2 use a linear parametrisation
$\bm{\theta}\mapsto\mathbf{K}_{\bm{\theta}}=\sum_{m=1}^{M}\theta_{m}\mathbf{K}_{m}$,
much of the results in MLM 2.3 and all the aforementioned discussion are still
applicable to non-linear parametrisations.

###  2.1 Multiple Kernel Learning

A widely adopted regularised risk principle, known as _multiple kernel
learning_ (MKL) (Christianini et al., 2001; Lanckriet et al., 2004; Bach et
al., 2004), is to minimise equation (1) w.r.t. the kernel parameters
$\bm{\theta}$ as well. One obvious caveat is that for any fixed $\mathbf{u}$,
equation (1) becomes ever smaller as $\theta_{m}\to\infty$: it cannot per se
play a meaningful statistical role. In order to prevent this, researchers
constrain the domain of $\bm{\theta}\in\bm{\Theta}$ and obtain

$\min_{\bm{\theta}\in\bm{\Theta}}\min_{\mathbf{u}}\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}+\ell(\mathbf{y},\mathbf{u}),$

where $\bm{\Theta}=\\{\bm{\theta}\succeq\mathbf{0},\>\|\bm{\theta}\|_{2}\leq
1\\}$ or
$\bm{\Theta}=\\{\bm{\theta}\succeq\mathbf{0},\>\mathbf{1}^{\top}\bm{\theta}\leq
1\\}$ (Varma and Ray, 2007). Notably, these constraints are imposed
independently of the statistical problem, the model and of the parametrization
$\bm{\theta}\mapsto\mathbf{K}_{\bm{\theta}}$. The Lagrangian form of the MKL
problem with parameter $\lambda$ and a general $p$-norm unit ball constraint
where $p\geq 1$ (Kloft et al., 2009) is given by

$\min_{\bm{\theta}\succeq\mathbf{0}}\phi_{\text{MKL}}(\bm{\theta}),\;\text{
where
}\phi_{\text{MKL}}(\bm{\theta}):=\min_{\mathbf{u}}\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}+\ell(\mathbf{y},\mathbf{u})+\underbrace{\lambda\cdot\mathbf{1}^{\top}\bm{\theta}^{p}}_{\rho(\bm{\theta})},\;\lambda>0.$
(3)

Since, the _regulariser_ $\rho(\bm{\theta})$ for the kernel parameter
$\bm{\theta}$ is convex, the map
$(\mathbf{u},\mathbf{K})\mapsto\mathbf{u}^{\top}\mathbf{K}^{-1}\mathbf{u}$ is
jointly convex for $\mathbf{K}\succeq\mathbf{0}$ (Boyd and Vandenberghe, 2002)
and the parametrisation $\bm{\theta}\mapsto\mathbf{K}_{\bm{\theta}}$ is
linear, MKL is a jointly convex problem for $\bm{\theta}\succeq\mathbf{0}$
whenever the loss function $\ell(\mathbf{y},\mathbf{u})$ is convex.
Furthermore, there are efficient algorithms to solve equation (3) for large
models (Sonnenburg et al., 2006).

###  2.2 Joint MAP Estimation

Adopting a probabilistic MAP viewpoint, we can minimise equation (2) w.r.t.
$\mathbf{u}$ and $\bm{\theta}\succeq\mathbf{0}$:

$\min_{\bm{\theta}\succeq\mathbf{0}}\phi_{\text{MAP}}(\bm{\theta}),\;\text{
where
}\phi_{\text{MAP}}(\bm{\theta}):=\min_{\mathbf{u}}\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}-2\ln\mathbb{P}(\mathbf{y}|\mathbf{u})+\ln|\mathbf{K}_{\bm{\theta}}|.$
(4)

While equation (3) and equation (4) share the “inner solution”
$\hat{\mathbf{u}}$ for fixed $\mathbf{K}_{\bm{\theta}}$ – in case the loss
$\ell(\mathbf{y},\mathbf{u})$ corresponds to a likelihood
$\mathbb{P}(\mathbf{y}|\mathbf{u})$ – they are different when it comes to
optimising $\bm{\theta}$. The _joint MAP_ problem is not in general jointly
convex in $(\bm{\theta},\mathbf{u})$, since
$\bm{\theta}\mapsto\ln|\mathbf{K}_{\bm{\theta}}|$ is concave, see figure 2.
However, it is always a well-posed statistical procedure, since
$\ln|\mathbf{K}_{\bm{\theta}}|\to\infty$ as $\theta_{m}\to\infty$ for all $m$.

Figure 1: Convex upper bounds on (the concave non-decreasing)
$\ln|\mathbf{K}_{\bm{\theta}}|$

By Fenchel duality, we can represent any concave non-decreasing function and
hence the log determinant function by
$\ln|\mathbf{K}_{\bm{\theta}}|=\min_{\bm{\lambda}\succeq\mathbf{0}}\bm{\lambda}^{\top}|\bm{\theta}|^{p}-g^{*}(\bm{\lambda})$.
As a consequence, we obtain a piecewise polynomial upper bound for any
particular value $\bm{\lambda}$.

We show in the following, how the regularisers
$\rho(\bm{\theta})=\lambda\left\|\bm{\theta}\right\|_{p}^{p}$ of equation (3)
can be related to the probabilistic term
$f(\bm{\theta})=\ln|\mathbf{K}_{\bm{\theta}}|$. In fact, the same reasoning
can be applied to any concave non-decreasing function.

Since the function $\bm{\theta}\mapsto
f(\bm{\theta})=\ln|\mathbf{K}_{\bm{\theta}}|$, $\bm{\theta}\succeq\mathbf{0}$
is jointly concave, we can represent it by
$f(\bm{\theta})=\min_{\bm{\lambda}\succeq\mathbf{0}}\bm{\lambda}^{\top}\bm{\theta}-f^{*}(\bm{\lambda})$
where $f^{*}(\bm{\lambda})$ denotes Fenchel dual of $f(\bm{\theta})$.
Furthermore, the mapping
$\bm{\vartheta}\mapsto\ln|\sum_{m=1}^{M}\sqrt[p]{\vartheta_{m}}\mathbf{K}_{m}|=f(\sqrt[p]{\bm{\vartheta}})=g(\bm{\vartheta})$,
$\bm{\vartheta}\succeq\mathbf{0}$ is jointly concave due to the composition
rule (Boyd and Vandenberghe, 2002, §3.2.4), because
$\bm{\vartheta}\mapsto\sqrt[p]{\bm{\vartheta}}$ is jointly concave and
$\bm{\theta}\mapsto f(\bm{\theta})$ is non-decreasing in all components
$\theta_{m}$ as all matrices $\mathbf{K}_{m}$ are positive (semi-)definite
which guarantees that the eigenvalues of $\mathbf{K}_{\bm{\theta}}$ increase
as $\theta_{m}$ increases. Thus we can – similarly to Zhang (2010) – represent
$\ln|\mathbf{K}_{\bm{\theta}}|$ as

$\ln|\mathbf{K}_{\bm{\theta}}|=f(\bm{\theta})=g(\bm{\vartheta})=\min_{\bm{\lambda}\succeq\mathbf{0}}\bm{\lambda}^{\top}\bm{\vartheta}-g^{*}(\bm{\lambda})=\min_{\bm{\lambda}\succeq\mathbf{0}}\bm{\lambda}^{\top}|\bm{\theta}|^{p}-g^{*}(\bm{\lambda}).$

Choosing a particular value $\bm{\lambda}=\lambda\cdot\mathbf{1}$, we obtain
the bound
$\ln|\mathbf{K}_{\bm{\theta}}|\leq\lambda\cdot\left\|\bm{\theta}\right\|_{p}^{p}-g^{*}(\lambda\cdot\mathbf{1})$.
Figure 1 illustrates the bounds for $p=1$ and $p=2$. The bottom line is that
one can interpret the regularisers
$\rho(\bm{\theta})=\lambda\left\|\bm{\theta}\right\|_{p}^{p}$ in equation (3)
as corresponding to parametrised upper bounds to the
$\ln|\mathbf{K}_{\bm{\theta}}|$ part in equation (4), hence
$\phi_{\text{MKL}}(\bm{\theta})=\psi_{\text{MAP}}(\bm{\theta},\bm{\lambda}=\lambda\cdot\mathbf{1})$,
where
$\phi_{\text{MAP}}(\bm{\theta})=\min_{\bm{\lambda}\succeq\mathbf{0}}\psi_{\text{MAP}}(\bm{\theta},\bm{\lambda})$.
Far from an ad hoc choice to keep $\bm{\theta}$ small, the
$\ln|\mathbf{K}_{\bm{\theta}}|$ term embodies the Occam’s razor concept behind
MAP estimation: overly large $\bm{\theta}$ are ruled out, since their
explanation of the data $\mathbf{y}$ is extremely unlikely under the prior
$\mathbb{P}(\mathbf{u};\bm{\theta})$. The Occam’s razor effect depends
crucially on the proper normalization of the prior (MacKay, 1992). For
example, the weighting parameter $C$ of $k$ ($k=C\tilde{k}$) can be learned by
joint MAP: if $C=e^{c}$, then equation (4) is convex in $c$ for any fixed
$\mathbf{u}$. If kernel-regularised estimation equation (1) is interpreted as
MAP estimation under a GP prior equation (2), the correct extension to kernel
learning is joint MAP: the MKL criterion equation (3) lacks prior
normalization, which renders MAP w.r.t. $\bm{\theta}$ meaningful in the first
place. From a non-probabilistic viewpoint, the $\ln|\mathbf{K}_{\bm{\theta}}|$
term comes with a model and data dependent structure at least as complex as
the rest of equation (3).

While the MKL objective, equation (3), enjoys the benefit of being convex in
the (linear) kernel parameters $\bm{\theta}$, this does not hold true for
joint MAP estimation, equation (4), in general. We illustrate the differences
in figure 2. The function $\psi_{\text{MAP}}(\bm{\theta},\mathbf{u})$ is a
building block of the MAP objective
$\phi_{\text{MAP}}(\bm{\theta})=\min_{\mathbf{u}}[\psi_{\text{MAP}}(\bm{\theta},\mathbf{u})-2\ln\mathbb{P}(\mathbf{y}|\mathbf{u})]$,
where

$\psi_{\text{MAP}}(\bm{\theta},\mathbf{u})=\underbrace{\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}}_{\psi_{\cup}(\bm{\theta},\mathbf{u})}+\underbrace{\ln|\mathbf{K}_{\bm{\theta}}|}_{\psi_{\cap}(\bm{\theta})}\leq\psi_{\text{MKL}}(\bm{\theta},\mathbf{u})-g^{*}(\lambda\cdot\mathbf{1}),\>\psi_{\text{MKL}}(\bm{\theta},\mathbf{u})=\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}+\lambda\left\|\bm{\theta}\right\|_{p}^{p}.$

More concretely, $\psi_{\text{MAP}}(\bm{\theta},\mathbf{u})$ is a sum of a
nonnegative, jointly convex function $\psi_{\cup}(\bm{\theta},\mathbf{u})$
that is strictly decreasing in every component $\theta_{m}$ and a concave
function $\psi_{\cap}(\bm{\theta})$ that is strictly increasing in every
component $\theta_{m}$. Both functions $\psi_{\cup}(\bm{\theta},\mathbf{u})$
and $\psi_{\cap}(\bm{\theta})$ alone do not have a stationary point due to
their monotonicity in $\theta_{m}$. However, their sum can have (even
multiple) stationary points as shown in figure 2 on the left left. We can
show, that the map
$\mathbf{K}\mapsto\mathbf{u}^{\top}\mathbf{K}^{-1}\mathbf{u}+\ln|\mathbf{K}|$
is _invex_ i.e. every stationary point $\hat{\mathbf{K}}$ is a global minimum.
Using the convexity of
$\mathbf{A}\mapsto\mathbf{u}^{\top}\mathbf{A}\mathbf{u}-\ln|\mathbf{A}|$ (Boyd
and Vandenberghe, 2002) and the fact that the derivative of
$\mathbf{A}\mapsto\mathbf{A}^{-1}$ for $\mathbf{A}\in\mathbb{R}^{n\times n}$,
$\mathbf{A}\succ\mathbf{0}$ has full rank $n^{2}$, we see by Mishra and Giorgi
(2008, theorem 2.1) that
$\mathbf{K}\mapsto\mathbf{u}^{\top}\mathbf{K}^{-1}\mathbf{u}+\ln|\mathbf{K}|$
is indeed invex.

Often, the MKL objective for the case $p=1$ is motivated by the fact that the
optimal solution $\bm{\theta}^{\star}$ is _sparse_ (e.g. Sonnenburg et al.,
2006), meaning that many components $\theta_{m}$ are zero. Figure 2
illustrates that $\phi_{\text{MAP}}(\bm{\theta})$ also yields sparse
solutions; in fact it enforces even more sparsity. In MKL,
$\phi_{\text{MAP}}(\bm{\theta})$ is simply relaxed to a convex objective
$\phi_{\text{MKL}}(\bm{\theta})$ at the expense of having only a single less
sparse solution.

Figure 2: Convex and non-convex building blocks of the MKL and MAP objective
function

#### Intuition for the Gaussian Case

We can gain further intuition about the criteria $\phi_{\text{MKL}}$ and
$\phi_{\text{MAP}}$ by asking which _matrices_ $\mathbf{K}$ minimise them. For
simplicity, assume that
$\mathbb{P}(\mathbf{y}|\mathbf{u})=\mathcal{N}(\mathbf{y}|\mathbf{u},\sigma^{2}\mathbf{I})$
and $n/C=\sigma^{2}$, hence
$\ell(\mathbf{y},\mathbf{u})=\frac{1}{\sigma^{2}}\left\|\mathbf{y}-\mathbf{u}\right\|_{2}^{2}$.
The inner minimiser $\hat{\mathbf{u}}$ for both $\phi_{\text{MKL}}$ and
$\phi_{\text{MAP}}$ is given by
$\mathbf{K}_{\bm{\theta}}^{-1}\hat{\mathbf{u}}=(\mathbf{K}_{\bm{\theta}}+\sigma^{2}\mathbf{I})^{-1}\mathbf{y}$.
With $\sigma^{2}\to 0$, we find for joint MAP that
$\frac{\partial}{\partial\mathbf{K}}\phi_{\text{MAP}}=\mathbf{0}$ results in
$\hat{\mathbf{K}}=\mathbf{y}\mathbf{y}^{\top}$. While this “nonparametric”
estimate requires smoothing to be useful in practice, closeness to
$\mathbf{y}\mathbf{y}^{\top}$ is fundamental to covariance estimation and can
be found in regularised risk kernel learning work (Christianini et al., 2001).
On the other hand, for $\text{tr}(\mathbf{K}_{m})=1$ and hence
$\rho(\bm{\theta})=\lambda\text{tr}(\mathbf{K}_{\bm{\theta}})=\lambda\left\|\bm{\theta}\right\|_{1}$,
$\frac{\partial}{\partial\mathbf{K}}\phi_{\text{MKL}}=\mathbf{0}$ leads to
$\hat{\mathbf{K}}^{2}=\lambda^{-1}\mathbf{y}\mathbf{y}^{\top}$: an odd way of
estimating covariance, not supported by any statistical literature we are
aware of.

###  2.3 Marginal Likelihood Maximisation

While the joint MAP criterion uses a properly normalised prior distribution,
it is still not probabilistically consistent. Kernel learning amounts to
finding a value $\hat{\bm{\theta}}$ of high data likelihood, no matter what
the latent function $u(\cdot)$ is. The correct likelihood to be maximised is
_marginal_ :
$\mathbb{P}(\mathbf{y}|\bm{\theta})=\int\mathbb{P}(\mathbf{y}|\mathbf{u})\mathbb{P}(\mathbf{u}|\bm{\theta})\text{d}\mathbf{u}$
(“max-sum”), while joint MAP employs the plug-in surrogate
$\max_{\mathbf{u}}\mathbb{P}(\mathbf{y}|\mathbf{u})\mathbb{P}(\mathbf{u}|\bm{\theta})$
(“max-max”). _Marginal likelihood maximisation_ (MLM) is also known as
Bayesian estimation, and it underlies the EM algorithm or maximum likelihood
learning of conditional random fields just as well: complexity is controlled
(and overfitting avoided) by averaging over unobserved variables $\mathbf{u}$
(MacKay, 1992), rather than plugging in some point estimate $\hat{\mathbf{u}}$

$\phi_{\text{MLM}}(\bm{\theta}):=-2\ln\int\mathcal{N}(\mathbf{u}|\mathbf{0},\mathbf{K}_{\bm{\theta}})\mathbb{P}(\mathbf{y}|\mathbf{u})\text{d}\mathbf{u}.$
(5)

#### The Gaussian Case

Before developing a general MLM approximation, we note an important
analytically solvable exception: for Gaussian likelihood
$\mathbb{P}(\mathbf{y}|\mathbf{u})=\mathcal{N}(\mathbf{y}|\mathbf{u},\sigma^{2}\mathbf{I})$,
$\mathbb{P}(\mathbf{y}|\bm{\theta})=\mathcal{N}(\mathbf{y}|\mathbf{0},\mathbf{K}_{\bm{\theta}}+\sigma^{2}\mathbf{I})$,
and MLM becomes

$\phi_{\text{GAU}}(\bm{\theta}):=\mathbf{y}^{\top}(\mathbf{K}_{\bm{\theta}}+\sigma^{2}\mathbf{I})^{-1}\mathbf{y}+\ln|\mathbf{K}_{\bm{\theta}}+\sigma^{2}\mathbf{I}|.$
(6)

Even if the primary purpose is classification, the Gaussian likelihood is used
for its analytical simplicity (Kapoor et al., 2009). Only for the Gaussian
case, joint MAP and MLM have an analytically closed form. From the product
formula of Gaussians (Brookes, 2005, §5.1)

$\mathbb{Q}(\mathbf{u}):=\mathcal{N}(\mathbf{u}|\mathbf{0},\mathbf{K}_{\bm{\theta}})\mathcal{N}(\mathbf{y}|\mathbf{u},\bm{\Gamma})=\mathcal{N}(\mathbf{y}|\mathbf{0},\mathbf{K}_{\bm{\theta}}+\bm{\Gamma})\mathcal{N}(\mathbf{u}|\mathbf{m},\mathbf{V}),$

where $\mathbf{V}=(\mathbf{K}_{\bm{\theta}}^{-1}+\bm{\Gamma}^{-1})^{-1}$ and
$\mathbf{m}=\mathbf{V}\bm{\Gamma}^{-1}\mathbf{y}$ we can deduce that

$-2\ln\int\mathbb{Q}(\mathbf{u})\text{d}\mathbf{u}=\ln|\mathbf{K}_{\bm{\theta}}^{-1}+\bm{\Gamma}^{-1}|+\min_{\mathbf{u}}[-2\ln\mathbb{Q}(\mathbf{u})]-n\ln|2\pi|.$
(7)

Using $\sigma^{2}\mathbf{I}=\bm{\Gamma}$ and
$\min_{\mathbf{u}}[-2\ln\mathbb{Q}(\mathbf{u})]=-2\ln\mathbb{Q}(\mathbf{m})$,
we see that by

$\phi_{\text{MAP/GAU}}(\bm{\theta}):\stackrel{{\scriptstyle\text{c}}}{{=}}\phi_{\text{GAU}}(\bm{\theta})-\ln|\mathbf{K}_{\bm{\theta}}^{-1}+\sigma^{-2}\mathbf{I}|\stackrel{{\scriptstyle\text{c}}}{{=}}\mathbf{y}^{\top}(\mathbf{K}_{\bm{\theta}}+\sigma^{2}\mathbf{I})^{-1}\mathbf{y}+\ln|\mathbf{K}_{\bm{\theta}}|$
(8)

MLM and MAP are very similar for the Gaussian case.

The “ridge regression” approximation is also used together with $p$-norm
constraints instead of the $\ln|\mathbf{K}_{\bm{\theta}}|$ term (Cortes et
al., 2009)

$\phi_{\text{RR}}(\bm{\theta}):=\mathbf{y}^{\top}(\mathbf{K}_{\bm{\theta}}+\sigma^{2}\mathbf{I})^{-1}\mathbf{y}+\lambda\left\|\bm{\theta}\right\|_{p}^{p}.$
(9)

Unfortunately, most GP methods to date work with a Gaussian likelihood for
simplicity, a restriction which often proves short-sighted. Gaussian-linear
models come with unrealistic properties, and benefits of MLM over joint MAP
cannot be realised.

Kernel parameter learning has been an integral part of probabilistic GP
methods from the very beginning. Williams and Rasmussen (1996) proposed MLM
for Gaussian noise equation 6, fifteen years ago. They treated sums of
exponential and linear kernels as well as learning lengthscales (ARD),
predating recent proposals such as “products of kernels” (Varma and Babu,
2009).

#### The General Case

In general, joint MAP always has the analytical form equation 4, while
$\mathbb{P}(\mathbf{y}|\bm{\text{$\theta$}})$ can only be approximated. For
non-Gaussian $\mathbb{P}(\mathbf{y}|\mathbf{u})$, numerous approximate
inference methods have been proposed, specifically motivated by learning
kernel parameters via MLM. The simplest such method is Laplace’s
approximation, applied to GP binary and multi-way classification by Williams
and Barber (1998): starting with convex joint MAP,
$\ln\mathbb{P}(\mathbf{y},\mathbf{u})$ is expanded to second order around the
posterior mode $\hat{\mathbf{u}}$. More recent approximations Girolami and
Rogers (2005); Girolami and Zhong (2006) can be much more accurate, yet come
with non-convex problems and less robust algorithms (Nickisch and Rasmussen,
2008). In this paper, we concentrate on the variational lower bound relaxation
(VB) by Jaakkola and Jordan (2000), which is convex for log-concave
likelihoods $\mathbb{P}(\mathbf{y}|\mathbf{u})$ (Nickisch and Seeger, 2009),
providing a novel simple and efficient algorithm. While our VB approximation
to MLM is more expensive to run than joint MAP for non-Gaussian likelihood
(even using Laplace’s approximation), the implementation complexity of our VB
algorithm is comparable to what is required in the Gaussian noise case
equation 6.

More, specifically, we exploit that super-Gaussian of likelihoods
$\mathbb{P}(y_{i}|u_{i})$ can be lower bounded by scaled Gaussians
$\mathcal{N}_{\gamma_{i}}$ of any width $\gamma_{i}$:

$\mathbb{P}(y_{i}|u_{i})=\max_{\gamma_{i}>0}\mathcal{N}_{\gamma_{i}}=\max_{\gamma_{i}>0}\exp\left(\beta_{i}u_{i}-\frac{u_{i}^{2}}{2\gamma_{i}}-\frac{1}{2}h_{i}(\gamma_{i})\right),$

where $\beta_{i}\propto y_{i}$ are constants, and $h_{i}(\cdot)$ is convex
(Nickisch and Seeger, 2009) whenever the likelihood $\mathbb{P}(y_{i}|u_{i})$
is log-concave. If the posterior distribution is
$\mathbb{P}(\mathbf{u}|\mathbf{y})=Z^{-1}\mathbb{P}(\mathbf{y}|\mathbf{u})\mathbb{P}(\mathbf{u})$,
then $\ln Z\geq Ce^{-\psi_{\text{VB}}(\bm{\theta},\bm{\gamma})/2}$ by plugging
in these bounds, where $C$ is a constant and

$\phi_{\text{VB}}(\bm{\theta}):=\min_{\bm{\gamma}\succ\mathbf{0}}\psi_{\text{VB}}(\bm{\theta},\bm{\gamma}),\quad\psi_{\text{VB}}(\bm{\theta},\bm{\gamma}):=h(\bm{\gamma})-2\ln\int\mathcal{N}(\mathbf{u}|\mathbf{0},\mathbf{K}_{\bm{\theta}})e^{\mathbf{u}^{\top}(\bm{\beta}-\frac{1}{2}\bm{\Gamma}{}^{-1}\mathbf{u})}\text{d}\mathbf{u},$
(10)

$h(\bm{\gamma}):=\sum_{i}h_{i}(\gamma_{i})$,
$\bm{\Gamma}:=\text{dg}(\bm{\gamma})$. The variational
relaxation111Generalisations to other super-Gaussian potentials (log-concave
or not) or models including linear couplings and mixed potentials are given by
Nickisch and Seeger (2009). amounts to maximising the lower bound, which means
that $\mathbb{P}(\mathbf{u}|\mathbf{y})$ is fitted by the _Gaussian_
approximation $\mathbb{Q}(\mathbf{u}|\mathbf{y};\bm{\gamma})$ with covariance
matrix $\mathbf{V}=(\mathbf{K}_{\bm{\theta}}^{-1}+\bm{\Gamma}{}^{-1})^{-1}$
(Nickisch and Seeger, 2009). Alternatively, we can interpret
$\psi_{\text{VB}}(\bm{\theta},\bm{\gamma})$ as an upper bound on the Kullback-
Leibler divergence
$\text{KL}(\mathbb{Q}(\mathbf{u}|\mathbf{y};\bm{\gamma})||\mathbb{P}(\mathbf{u}|\mathbf{y}))$
(Nickisch, 2010, §2.5.9), a measure for the dissimilarity between the exact
posterior $\mathbb{P}(\mathbf{u}|\mathbf{y})$ and the parametrised Gaussian
approximation $\mathbb{Q}(\mathbf{u}|\mathbf{y};\bm{\gamma})$.

Finally, note that by equation (7),
$\psi_{\text{VB}}(\bm{\theta},\bm{\gamma})$ can also be written as

$\psi_{\text{VB}}(\bm{\theta},\bm{\gamma})=\ln|\mathbf{K}_{\bm{\theta}}^{-1}+\bm{\Gamma}^{-1}|+h(\bm{\gamma})+\min_{\mathbf{u}}R(\mathbf{u},\bm{\theta},\bm{\gamma})+\ln|\mathbf{K}_{\bm{\theta}}|,$
(11)

where
$R(\mathbf{u},\bm{\theta},\bm{\gamma})=\mathbf{u}^{\top}(\mathbf{K}_{\bm{\theta}}^{-1}+\bm{\Gamma}^{-1})\mathbf{u}-2\bm{\beta}^{\top}\mathbf{u}$.
Using the concavity of
$\bm{\gamma}^{-1}\mapsto\ln|\mathbf{K}_{\bm{\theta}}^{-1}+\bm{\Gamma}^{-1}|$
and Fenchel duality
$\ln|\mathbf{K}_{\bm{\theta}}^{-1}+\bm{\Gamma}^{-1}|=\min_{\mathbf{z}\succ\mathbf{0}}\mathbf{z}^{\top}\bm{\gamma}^{-1}-g_{\bm{\theta}}^{*}(\mathbf{z})=\hat{\mathbf{z}}_{\bm{\theta}}^{\top}\bm{\gamma}^{-1}-g_{\bm{\theta}}^{*}(\hat{\mathbf{z}}_{\bm{\theta}})$,
with the optimal value $\hat{\mathbf{z}}_{\bm{\theta}}=\text{dg}(\mathbf{V})$,
we can reformulate $\psi_{\text{VB}}(\bm{\theta},\bm{\gamma})$ as

$\psi_{\text{VB}}(\bm{\theta},\bm{\gamma})=\min_{\mathbf{z}\succ\mathbf{0}}[\mathbf{z}^{\top}\bm{\gamma}^{-1}-g_{\bm{\theta}}^{*}(\mathbf{z})]+h(\bm{\gamma})+\min_{\mathbf{u}}R(\mathbf{u},\bm{\theta},\bm{\gamma})+\ln|\mathbf{K}_{\bm{\theta}}|,$

which allows to perform the minimisation w.r.t. $\bm{\gamma}$ in closed form
(Nickisch, 2010, §3.5.6):

$\phi_{\text{VB}}(\bm{\theta})=\min_{\mathbf{z}\succ\mathbf{0}}\psi_{\text{VB}}(\bm{\theta},\mathbf{z}),\quad\psi_{\text{VB}}(\bm{\theta},\mathbf{z})=\min_{\mathbf{u}}\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}+\tilde{\ell}_{\mathbf{z}}(\mathbf{y},\mathbf{u})-g_{\bm{\theta}}^{*}(\mathbf{z})+\ln|\mathbf{K}_{\bm{\theta}}|,$
(12)

where
$\tilde{\ell}_{\mathbf{z}}(\mathbf{y},\mathbf{u}):=2\bm{\beta}^{\top}(\mathbf{v}-\mathbf{u})-2\ln\mathbb{P}(\mathbf{y}|\mathbf{v})$
and finally
$\mathbf{v}=\text{sign}(\mathbf{u})\odot\sqrt{\mathbf{u}^{2}+\mathbf{z}}$.
Note that for $\mathbf{z}=\mathbf{0}$, we exactly recover joint MAP
estimation, equation (4), as $\mathbf{z}=\mathbf{0}$ implies
$\mathbf{u}=\mathbf{v}$ and
$\tilde{\ell}_{\mathbf{z}}(\mathbf{y},\mathbf{u})=\ell(\mathbf{y},\mathbf{u})$.
For fixed $\bm{\theta}$, the optimal value
$\hat{\mathbf{z}}_{\bm{\theta}}=\text{dg}(\mathbf{V})$ corresponds to the
marginal variances of the Gaussian approximation
$\mathbb{Q}(\mathbf{u}|\mathbf{y};\bm{\gamma})$: Variational inference
corresponds to variance-smoothed joint MAP estimation (Nickisch, 2010) with a
loss function $\tilde{\ell}(\mathbf{y},\mathbf{u},\bm{\theta})$ that
explicitly depends on the kernel parameters $\bm{\theta}$. We have two
equivalent representations of the loss
$\tilde{\ell}(\mathbf{y},\mathbf{u},\bm{\theta})$ that directly follow from
equations (11) and (12):

$\displaystyle\tilde{\ell}(\mathbf{y},\mathbf{u},\bm{\theta})$
$\displaystyle=$
$\displaystyle\min_{\bm{\gamma}\succ\mathbf{0}}[\ln|\mathbf{K}_{\bm{\theta}}^{-1}+\bm{\Gamma}^{-1}|+h(\bm{\gamma})+\mathbf{u}^{\top}\bm{\Gamma}^{-1}\mathbf{u}-2\bm{\beta}^{\top}\mathbf{u}],\>\text{and}$
$\displaystyle\tilde{\ell}(\mathbf{y},\mathbf{u},\bm{\theta})$
$\displaystyle=$
$\displaystyle\min_{\mathbf{z}\succ\mathbf{0}}[2\bm{\beta}^{\top}(\mathbf{v}-\mathbf{u})-2\ln\mathbb{P}(\mathbf{y}|\mathbf{v})-g_{\bm{\theta}}^{*}(\mathbf{z})],\;\mathbf{v}=\text{sign}(\mathbf{u})\odot\sqrt{\mathbf{u}^{2}+\mathbf{z}}.$

Our VB problem is
$\min_{\bm{\theta}\succeq\mathbf{0},\bm{\gamma}\succ\mathbf{0}}\psi_{\text{VB}}(\bm{\theta},\bm{\gamma})$
or equivalently
$\min_{\bm{\theta}\succeq\mathbf{0},\mathbf{z}\succ\mathbf{0}}\psi_{\text{VB}}(\bm{\theta},\mathbf{z})$.
The inner variables here are $\bm{\gamma}$ and $\mathbf{z}$, in addition to
$\mathbf{u}$ in joint MAP. There are further similarities: since
$\psi_{\text{VB}}(\bm{\theta},\bm{\gamma})=-2\ln\int
e^{-R(\mathbf{u},\bm{\gamma},\bm{\theta})}\text{d}\mathbf{u}+h(\bm{\gamma})+\ln|2\pi\mathbf{K}_{\bm{\theta}}|$,
$(\bm{\gamma},\bm{\theta})\mapsto\psi_{\text{VB}}-\ln|\mathbf{K}_{\bm{\theta}}|$
is jointly convex for $\bm{\gamma}\succ\mathbf{0}$,
$\bm{\theta}\succeq\mathbf{0}$, by the joint convexity of
$(\mathbf{u},\bm{\gamma},\bm{\theta})\mapsto R$ and Prékopa’s theorem (Boyd
and Vandenberghe, 2002, §3.5.2). Joint MAP and VB share the same convexity
structure. In contrast, approximating $\mathbb{P}(\mathbf{y}|\bm{\theta})$ by
other techniques like Expectation Propagation (Minka, 2001) or general
Variational Bayes (Opper and Archambeau, 2009) does not even constitute convex
problems for fixed $\bm{\theta}$.

###  2.4 Summary and Taxonomy

Name | Objective function
---|---
Marginal Likelihood Maximisation | $\phi_{\text{MLM}}(\bm{\theta})=-2\ln\left[\int\mathcal{N}(\mathbf{u}|\mathbf{0},\mathbf{K}_{\bm{\theta}})\mathbb{P}(\mathbf{y}|\mathbf{u})\text{d}\mathbf{u}\right]$
Variational Bounds |  $\phi_{\text{VB}}(\bm{\theta})=\min_{\bm{\gamma}\succ\mathbf{0}}\psi_{\text{VB}}(\bm{\theta},\bm{\gamma})\geq\phi_{\text{MLM}}(\bm{\theta})$ by $\mathbb{P}(y_{i}|u_{i})\geq\mathcal{N}_{\gamma_{i}}$
Maximum A Posteriori | $\phi_{\text{MAP}}(\bm{\theta})=-2\ln\left[\max_{\mathbf{u}}\mathcal{N}(\mathbf{u}|\mathbf{0},\mathbf{K}_{\bm{\theta}})\mathbb{P}(\mathbf{y}|\mathbf{u})\right]=\psi_{\text{VB}}(\bm{\theta},\mathbf{z}=\mathbf{0})$
Multiple Kernel Learning | $\phi_{\text{MKL}}(\bm{\theta})=\phi_{\text{MAP}}(\bm{\theta})+\lambda\left\|\bm{\theta}\right\|_{p}^{p}-\ln|\mathbf{K}_{\bm{\theta}}|=\psi_{\text{MAP}}(\bm{\theta},\bm{\lambda}=\lambda\cdot\mathbf{1})$

| General $\mathbb{P}(y_{i}|u_{i})$ |  | Gaussian $\mathbb{P}(y_{i}|u_{i})$ |
---|---|---|---|---
|  $\phi_{\text{MLM}}(\bm{\theta})$, eq. (5) | $\longrightarrow$ |  $\phi_{\text{GAU}}(\bm{\theta})$, eq. (6) |
Super-Gaussian Bounding |  $\downarrow$ |  |  $\downarrow$ | Bound is tight
|  $\phi_{\text{VB}}(\bm{\theta})$, eq. (10) | $\longrightarrow$ |  $\phi_{\text{GAU}}(\bm{\theta})$, eq. (6) |
Maximum instead of integral |  $\downarrow$ |  |  $\downarrow$ | Mode $\equiv$ mean
|  $\phi_{\text{MAP}}(\bm{\theta})$, eq. (4) | $\longrightarrow$ |  $\phi_{\text{MAP/GAU}}(\bm{\theta})$, eq. (8) |
Bound $\ln|\mathbf{K}_{\bm{\theta}}|\leq\lambda\left\|\bm{\theta}\right\|_{p}^{p}-g^{*}(\lambda\mathbf{1})$ |  $\downarrow$ |  |  $\downarrow$ |
|  $\phi_{\text{MKL}}(\bm{\theta})$, eq. (3) | $\longrightarrow$ |  $\phi_{\text{RR}}(\bm{\theta})$, eq. (9) |

Table 2: Taxonomy of kernel learning objective functions

The upper table visualises the relationship between several kernel learning
objective functions for arbitrary likelihood/loss functions: Marginal
likelihood maximisation (MLM) can be bounded by variational bounds (VB) and
maximum a posteriori estimation (MAP) is a special case
$\mathbf{z}=\mathbf{0}$ thereof. Finally multiple kernel learning (MKL) can be
understood as an upper bound to the MAP estimation objective
$\bm{\lambda}=\lambda\cdot\mathbf{1}$. The lower table complements the upper
table by also covering the analytically important Gaussian case.

In the last paragraphs, we have detailed how a variety of kernel learning
approaches can be obtained from Bayesian marginal likelihood maximisation in a
sequence of nested upper bounding steps. Table 2.4 nicely illustrates how many
kernel learning objectives are related to each other – either by upper bounds
or by Gaussianity assumptions. We can clearly see, that
$\phi_{\text{VB}}(\bm{\theta})$ – as an upper bound to the negative log
marginal likelihood – can be seen as the mother function. For a special case,
$\mathbf{z}=\mathbf{0}$, we obtain joint maximum a posteriori estimation,
where the loss functions does not depend on the kernel parameters. Going
further, a particular instance $\bm{\lambda}=\lambda\cdot\mathbf{1}$ yields
the widely use multiple kernel learning objective that becomes convex in the
kernel parameters $\bm{\theta}$. In the following, we will concentrate on the
optimisation and computational similarities between the approaches.

##  3 Algorithms

In this section, we derive a simple, provably convergent and efficient
algorithm for MKL, joint MAP and VB. We use the Lagrangian form of equation
(3) and $\ell(\mathbf{y},\mathbf{u}):=-2\ln\mathbb{P}(\mathbf{y}|\mathbf{u})$:

$\displaystyle\psi_{\text{MKL}}(\bm{\theta},\mathbf{u})$ $\displaystyle=$
$\displaystyle\mathbf{u}^{\top}\mathbf{K}^{-1}\mathbf{u}+\qquad\;\ell(\mathbf{y},\mathbf{u})\qquad\qquad\qquad\qquad\qquad\;\,+\lambda\cdot\mathbf{1}^{\top}\bm{\theta},\>\lambda>0,$
$\displaystyle\psi_{\text{MAP}}(\bm{\theta},\mathbf{u})$ $\displaystyle=$
$\displaystyle\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}+\qquad\;\ell(\mathbf{y},\mathbf{u})\qquad\qquad\qquad\qquad\qquad\;\,+\ln|\mathbf{K}_{\bm{\theta}}|,\quad\text{and}$
$\displaystyle\psi_{\text{VB}}(\bm{\theta},\mathbf{u})$ $\displaystyle=$
$\displaystyle\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}+\min_{\mathbf{z}\succ\mathbf{0}}\left[\ell(\mathbf{y},\mathbf{v})+2\bm{\beta}^{\top}(\mathbf{v}-\mathbf{u})-g_{\bm{\theta}}^{*}(\mathbf{z})\right]+\ln|\mathbf{K}_{\bm{\theta}}|,$
$\displaystyle\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad\text{where}\;\mathbf{v}=\text{sign}(\mathbf{u})\odot\sqrt{\mathbf{u}^{2}+\mathbf{z}}.$

Many previous algorithms use alternating minimization, which is easy to
implement but tends to converge slowly. Both $\phi_{\text{VB}}$ and
$\phi_{\text{MAP}}$ are jointly convex up to the concave
$\bm{\theta}\mapsto\ln|\mathbf{K}_{\bm{\theta}}|$ part. Since
$\ln|\mathbf{K}_{\bm{\theta}}|=\min_{\bm{\lambda}\succ\mathbf{0}}\bm{\lambda}^{\top}\bm{\theta}-f^{*}(\bm{\lambda})$
(Legendre duality, Boyd and Vandenberghe, 2002), joint MAP becomes
$\min_{\bm{\lambda}\succ\mathbf{0},\bm{\theta}\succeq\mathbf{0},\mathbf{u}}\phi_{\bm{\lambda}}(\bm{\theta},\mathbf{u})$
with
$\phi_{\bm{\lambda}}:=\mathbf{u}^{\top}\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{u}+\ell(\mathbf{y},\mathbf{u})+\bm{\lambda}^{\top}\bm{\theta}-f^{*}(\bm{\lambda})$
which is jointly convex in $(\bm{\theta},\mathbf{u})$. Algorithm 1 iterates
between refits of $\bm{\lambda}$ and joint Newton updates of
$(\bm{\theta},\mathbf{u})$.

0: Criterion
$\psi_{\\#}(\bm{\theta},\mathbf{u})=\tilde{\psi}_{\\#}(\bm{\theta},\mathbf{u})+\ln|\mathbf{K}_{\bm{\theta}}|$
to minimise for
$(\mathbf{u},\bm{\theta})\in\mathbb{R}^{n}\times\mathbb{R}_{+}^{M}$.

repeat

Newton $\min_{\mathbf{u}}\psi_{\\#}$ for fixed $\bm{\theta}$ (optional; few
steps).

Refit upper bound:
$\bm{\lambda}\leftarrow\nabla_{\bm{\theta}}\ln|\mathbf{K}_{\bm{\theta}}|=[\text{tr}(\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{K}_{1}),..,\text{tr}(\mathbf{K}_{\bm{\theta}}^{-1}\mathbf{K}_{M})]^{\top}$.

Compute joint Newton search direction $\mathbf{d}$ for
$\psi_{\bm{\lambda}}:=\tilde{\psi}_{\\#}+\bm{\lambda}^{\top}\bm{\theta}$:
$\nabla_{[\bm{\theta};\mathbf{u}]}^{2}\psi_{\bm{\lambda}}\mathbf{d}=-\nabla_{[\bm{\theta};\mathbf{u}]}\psi_{\bm{\lambda}}$.

Linesearch: Minimise $\psi_{\\#}(\alpha)$ i.e.
$\psi_{\\#}(\bm{\theta},\mathbf{u})$ along
$[\bm{\theta};\mathbf{u}]+\alpha\mathbf{d}$, $\alpha>0$.

until Outer loop converged

Algorithm 1 Double loop algorithm for joint MAP, MKL and VB.

The Newton direction costs $O(n^{3}+M\,n^{2})$, with $n$ the number of data
points and $M$ the number of base kernels. All algorithms discussed in this
paper require $O(n^{3})$ time, apart from the requirement to store the base
matrices $\mathbf{K}_{m}$. The convergence proof hinges on the fact that
$\phi$ and $\phi_{\bm{\lambda}}$ are tangentially equal (Nickisch and Seeger,
2009). Equivalently, the algorithm can be understood as Newton’s method, yet
dropping the part of the Hessian corresponding to the $\ln|\mathbf{K}|$ term
(note that
$\nabla_{(\mathbf{u},\bm{\theta})}\phi_{\bm{\lambda}}=\nabla_{(\mathbf{u},\bm{\theta})}\phi$
for the Newton direction computation). Exact Newton for MKL.

In practice, we use
$\mathbf{K}_{\bm{\theta}}=\sum_{m}\theta_{m}\mathbf{K}_{m}+\varepsilon\mathbf{I},\>\varepsilon=10^{-8}$
to avoid numerical problems when computing $\bm{\lambda}$ and
$\ln|\mathbf{K}_{\bm{\theta}}|$. We also have to enforce
$\bm{\theta}\succeq\mathbf{0}$ in algorithm 1, which is done by the barrier
method (Boyd and Vandenberghe, 2002). We minimise
$t\phi+\mathbf{1}^{\top}(\ln\bm{\theta})$ instead of $\phi$, increasing $t>0$
every few outer loop iterations.

A variant algorithm 1 can be used to solve VB in a different parametrisation
($\bm{\gamma}\succ\mathbf{0}$ replaces $\mathbf{u}$), which has the same
convexity structure as joint MAP. Transforming equation (10) similarly to
equation (6), we obtain

$\phi_{\text{VB}}(\bm{\theta})=\min_{\bm{\gamma}\succ\mathbf{0}}\ln|\mathbf{C}|-\ln|\bm{\Gamma}|+\bm{\beta}^{\top}\bm{\Gamma}\mathbf{C}^{-1}\bm{\Gamma}\bm{\beta}-\bm{\beta}^{\top}\bm{\Gamma}\bm{\beta}+h(\bm{\gamma})$
(13)

with $\mathbf{C}:=\mathbf{K}_{\bm{\theta}}+\bm{\Gamma}$, computed using the
Cholesky factorisation $\mathbf{C}=\mathbf{L}\mathbf{L}{}^{\top}$. They cost
$O(M\,n^{3})$ to compute, which is more expensive than for joint MAP or MKL.
Note that the cost $O(M\,n^{3})$ is not specific to our particular relaxation
or algorithm e.g. the Laplace MLM approximation (Williams and Barber, 1998),
solved using gradients w.r.t. $\bm{\theta}$ only, comes with the same
complexity.

##  4 Conclusion

We presented a unifying probabilistic viewpoint to multiple kernel learning
that derives regularised risk approaches as special cases of approximate
Bayesian inference. We provided an efficient and provably convergent
optimisation algorithm suitable for regression, robust regression and
classification.

Our taxonomy of multiple kernel learning approaches connected many previously
only loosely related ideas and provided insights into the common structure of
the respective optimisation problems. Finally, we proposed an algorithm to
solve the latter efficiently.

## References

  * Bach et al. (2004) Francis Bach, Gert Lanckriet, and Michael Jordan.  Multiple kernel learning, conic duality, and the SMO algorithm.  In _ICML_ , 2004.
  * Boyd and Vandenberghe (2002) Stephen Boyd and Lieven Vandenberghe.  _Convex Optimization_.  Cambridge University Press, 2002.
  * Brookes (2005) Mike Brookes.  The matrix reference manual, 2005.  URL http://www.ee.ic.ac.uk/hp/staff/dmb/matrix/intro.html.
  * Christianini et al. (2001) Nello Christianini, John Shawe-Taylor, André Elisseeff, and Jaz Kandola.  On kernel-target alignment.  In _NIPS_ , 2001.
  * Cortes et al. (2009) Corinna Cortes, Mehryar Mohri, and Afshin Rostamizadeh.  L2 regularization for learning kernels.  In _UAI_ , 2009.
  * Girolami and Rogers (2005) Mark Girolami and Simon Rogers.  Hierarchic Bayesian models for kernel learning.  In _ICML_ , 2005.
  * Girolami and Zhong (2006) Mark Girolami and Mingjun Zhong.  Data integration for classification problems employing Gaussian process.  In _NIPS_ , 2006.
  * Jaakkola and Jordan (2000) Tomi Jaakkola and Michael Jordan.  Bayesian parameter estimation via variational methods.  _Statistics and Computing_ , 10:25–37, 2000.
  * Kapoor et al. (2009) Ashish Kapoor, Kristen Grauman, Raquel Urtasun, and Trevor Darrell.  Gaussian processes for object categorization.  _IJCV_ , 2009.  doi: 10.1007/s11263-009-0268-3.
  * Kloft et al. (2009) Marius Kloft, Ulf Brefeld, Sören Sonnenburg, Pavel Laskov, Klaus-Robert Müller, and Alexander Zien.  Efficient and accurate lp-norm multiple kernel learning.  In _NIPS_ , 2009.
  * Lanckriet et al. (2004) Gert R. G. Lanckriet, Nello Cristianini, Peter Bartlett, Laurent El Ghaoui, and Michal I. Jordan.  Learning the kernel matrix with semidefinite programming.  _JMLR_ , 5:27–72, 2004.
  * MacKay (1992) Davic MacKay.  Bayesian interpolation.  _Neural Computation_ , 4(3):415–447, 1992.
  * Minka (2001) Tom Minka.  Expectation propagation for approximate Bayesian inference.  In _UAI_ , 2001.
  * Mishra and Giorgi (2008) Shashi Kant Mishra and Giorgio Giorgi.  _Invexity and optimization_.  Springer, 2008.
  * Nickisch (2010) Hannes Nickisch.  _Bayesian Inference and Experimental Design for Large Generalised Linear Models_.  PhD thesis, TU Berlin, 2010.
  * Nickisch and Rasmussen (2008) Hannes Nickisch and Carl Edward Rasmussen.  Approximations for binary Gaussian process classification.  _JMLR_ , 9:2035–2078, 2008.
  * Nickisch and Seeger (2009) Hannes Nickisch and Matthias Seeger.  Convex variational Bayesian inference for large scale generalized linear models.  In _ICML_ , 2009.
  * Opper and Archambeau (2009) Manfred Opper and Cédric Archambeau.  The variational Gaussian approximation revisited.  _Neural Computation_ , 21(3):786–792, 2009.
  * Rasmussen and Williams (2006) Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian Processes for Machine Learning_.  MIT Press, 2006.
  * Schölkopf and Smola (2002) Bernhard Schölkopf and Alex Smola.  _Learning with Kernels_.  MIT Press, 1st edition, 2002.
  * Sollich (2000) Peter Sollich.  Probabilistic methods for support vector machines.  In _NIPS_ , 2000.
  * Sonnenburg et al. (2006) Sören Sonnenburg, Gunnar Rätsch, Christin Schäfer, and Bernhard Schölkopf.  Large scale multiple kernel learning.  _JMLR_ , 7:1531–1565, 2006.
  * Varma and Babu (2009) Manik Varma and Bodla Rakesh Babu.  More generality in efficient multiple kernel learning.  In _ICML_ , 2009.
  * Varma and Ray (2007) Manik Varma and Debajyoti Ray.  Learning the discriminative power-invariance trade-off.  In _ICCV_ , 2007.
  * Williams and Barber (1998) Christopher K. I. Williams and David Barber.  Bayesian classification with Gaussian processes.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 20(12):1342–1351, 1998.
  * Williams and Rasmussen (1996) Christopher K. I. Williams and Carl Edward Rasmussen.  Gaussian processes for regression.  In _NIPS_ , 1996.
  * Zhang (2010) Tong Zhang.  Analysis of multi-stage convex relaxation for sparse regularization.  _JMLR_ , 11:1081–1107, 2010.