IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, IN PRESS

1

Identifying Corresponding Patches in SAR and
Optical Images with a Pseudo-Siamese CNN

arXiv:1801.08467v1 [eess.IV] 25 Jan 2018

Lloyd H. Hughes, Student Member, IEEE, Michael Schmitt, Senior Member, IEEE, Lichao Mou,
Yuanyuan Wang, Member, IEEE, and Xiao Xiang Zhu, Senior Member, IEEE

Abstract—This is the pre-print version, to read the final version
please go to IEEE Geoscience and Remote Sensing Letters on
IEEE Xplore.
In this letter, we propose a pseudo-siamese convolutional neural
network (CNN) architecture that enables to solve the task of
identifying corresponding patches in very-high-resolution (VHR)
optical and synthetic aperture radar (SAR) remote sensing
imagery. Using eight convolutional layers each in two parallel
network streams, a fully connected layer for the fusion of the
features learned in each stream, and a loss function based on
binary cross-entropy, we achieve a one-hot indication if two
patches correspond or not.
The network is trained and tested on an automatically generated dataset that is based on a deterministic alignment of
SAR and optical imagery via previously reconstructed and subsequently co-registered 3D point clouds. The satellite images, from
which the patches comprising our dataset are extracted, show
a complex urban scene containing many elevated objects (i.e.
buildings), thus providing one of the most difficult experimental
environments. The achieved results show that the network is
able to predict corresponding patches with high accuracy, thus
indicating great potential for further development towards a
generalized multi-sensor key-point matching procedure.
Index Terms—synthetic aperture radar (SAR), optical imagery,
data fusion, deep learning, convolutional neural networks (CNN),
image matching, deep matching

I. I NTRODUCTION

T

HE identification of corresponding image patches is
used extensively in computer vision and remote sensingrelated image analysis, especially in the framework of stereo
applications or co-registration issues. While many successful
hand-crafted approaches, specifically designed for the matching of optical images, exist [1], to this date the matching of
images acquired by different sensors still remains a widely
unsolved challenge [2]. This particularly holds for a joint
exploitation of SAR and optical imagery, caused by two
completely different sensing modalities: SAR imagery collects information about the physical properties of the scene
The authors are with Signal Processing in Earth Observation,
Technical University of Munich, 80333 Munich, Germany (email:
lloyd.hughes@tum.de,
m.schmitt@tum.de,
lichao.mou@dlr.de,
yuanyuan.wang@dlr.de, xiao.zhu@dlr.de)
X. Zhu is also with the Remote Sensing Technology Institute, German
Aerospace Center, 82234 Wessling, Germany
Manuscript received XXXX YY, ZZZZ.
This work was supported by the China Scholarship Council, the European
Research Council (ERC) under the European Unions Horizon 2020 research
and innovation programme (grant agreement No. ERC-2016-StG-714087,
Acronym: So2Sat), the Helmholtz Association under the framework of the
Young Investigators Group SiPEO (VH-NG-1018, www.sipeo.bgu.tum.de)
and the German Research Foundation (DFG), grant SCHM 3322/1-1.

and follows a range-based imaging geometry, while optical
imagery reflects the chemical characteristics of the scene and
follows a perspective imaging geometry. Hence, structures
elevated above the ground level, such as buildings or trees,
show strongly different appearances in both SAR and optical
images (cf. Fig. 1), in particular when dealing with very high
resolution data.

Fig. 1. A simple detached multi-story building as SAR amplitude image (left)
and optical photograph (right).

In order to deal with the problem of multi-sensor keypoint matching, several sophisticated approaches have been
proposed, e.g. exploiting phase congruency as a generalization
of gradient information [3]. However, even sophisticated handcrafted descriptors reach their limitations for highly resolving
data showing densely built-up urban scenes, which – in the
SAR case – is often difficult to interpret even for trained
experts.
Therefore, our work aims at learning a multi-sensor correspondence predictor for SAR and optical image patches
of state-of-the-art VHR data. Inspired by promising results
achieved in the context of stereo matching for optical imagery
[4], [5], we also make use of a convolutional neural network
(CNN). The major difference of our work to these purely
optical approaches is that we focus on the aforementioned,
distinctly more complicated multi-sensor setup, and therefore design a specific pseudo-siamese network architecture
with two separate, yet identical convolutional streams for
processing SAR and optical patches in parallel, instead of
a weight-shared siamese network in order to deal with the
heterogeneous nature of the input imagery.
II. T HE N ETWORK A RCHITECTURE
A. Pseudo-Siamese Convolutional Network
Since SAR and optical images lie on different manifolds,
it is not advisable to compare them directly by descriptors

IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, IN PRESS

IUT\E

2

IUT\E

IUT\E

IUT\E

IUT\E

UVZOIGRVGZIN

IUT\E

IUT\E

IUT\E

Details for Config.


^






^






 ^
^

 





^



^



^
















^



^






^


^



^






^ ^






IUT\





RGHKR









^






^





^ 
 ^





9'8VGZIN



Conv 2_2 + ReLU + BN

stride 2

Max Pooling 1

Max Pooling 1

conv. 3x3 (l2 norm), stride 1, padding 1

Conv 1_3 + ReLU + BN

Conv 2_3 + ReLU + BN

conv. 3x3 (l2 norm), stride 1, padding 1

Conv 1_4 + ReLU + BN

Conv 2_4 + ReLU + BN

stride 2

Max Pooling 2

Max Pooling 2

conv. 3x3 (l2 norm), stride 1, padding 1

Conv 1_5 + ReLU + BN

Conv 2_5 + ReLU + BN

conv. 3x3 (l2 norm), stride 1, padding 1

Conv 1_6 + ReLU + BN

Conv 2_6 + ReLU + BN

stride 2

Max Pooling 3

Max Pooling 3

conv. 3x3 (l2 norm), stride 1, padding 1

Conv 1_7 + ReLU + BN

Conv 2_7 + ReLU + BN

conv. 3x3 (l2 norm), stride 1, padding 1

Conv 1_8 + ReLU + BN

Conv 2_8 + ReLU + BN



conv. 3x3 (l2 norm), stride 2, padding 1

Conv 9 + ReLU + BN

conv. 3x3 (l2 norm), stride 2, padding 1

Conv 10 + ReLU + BN



^

IUT\E

Max Pooling 4
FC 1 + ReLU



IUT\E

Conv 2_1 + ReLU + BN

Conv 1_2 + ReLU + BN

stride 2



^
^
 



^

IUT\E




^

 ^




^



IUT\E






^




^


^









^




^ ^






IUT\E

Optical Conv. Stream

Conv 1_1 + ReLU + BN

conv. 3x3 (l2 norm), stride 1, padding 1

Concatenation Layer


^



^



^



^

IUT\E

LI

IUTIGZ



IUT\E

IUT\


^





^



^



^










LI

SAR Conv. Stream
conv. 3x3 (l2 norm), stride 1, padding 1

FC 2 + Softmax
Loss

IUT\E

Fig. 2. Proposed pseudo-Siamese network architecture and layer configuration.

designed for matching optical patches. Neither suitable are
conventional siamese CNN architectures (e.g. [5]), which are
designed to share their weights in each layer, as the input
data processed in the two feature extraction streams share
similar properties. In order to cope with the strongly different
geometric and radiometric appearances of SAR and optical
imagery, in [6] we proposed a pseudo-siamese network architecture with two separate, yet identical convolutional streams,
which process the SAR patch and the optical patch in parallel,
and only fuse the resulting information at a later decision
stage. Using this architecture, the network is constrained to
first learn meaningful representations of the input SAR patch
and the optical patch separately, and to combine them on a
higher level. The work presented in this letter is an extension
of [6] by improving the fusion part of the network architecture,
using a different training strategy, and resorting to non-locally
prefiltered SAR patches instead of temporal mean maps.
In addition, we evaluate the network on a deterministically
partitioned dataset instead of a randomly partitioned one, as
random partitioning will always cause positively biased results
due to overlapping regions in patches.
The architecture of the proposed network is shown in Fig. 2.
It is mainly inspired by the philosophy of the well-known
VGG Nets [7]. The SAR and optical image patches are passed
through a stack of convolutional layers, where we make use of
convolutional filters with a very small receptive field of 3 × 3,
rather than using larger ones, such as 5×5 or 7×7. The reason
is that 3×3 convolutional filters are the smallest kernels to capture patterns in different directions, such as center, up/down,
and left/right, but still have an advantage: the use of small
convolutional filters will increase the nonlinearities inside the
network and thus make the network more discriminative
The convolution stride in our network is fixed to 1 pixel;
the spatial padding of convolutional layer input is such that
the spatial resolution is preserved after convolution, i.e. the
padding is 1 pixel for the all 3 × 3 convolutional layers in our
network. Spatial pooling is achieved by carrying out seven
max-pooling layers which follow some of the convolutional
layers. They are used to reduce the dimensionality of the
feature maps. Max-pooling is performed over 2 × 2 pixel

windows with stride 2.
The fusion stage of our proposed network is made up of
two consecutive convolutional layers, followed by two fully
connected layers. The convolutional layers consist of 3 × 3
filters which operate over the concatenated feature maps of the
SAR and optical streams, in order to learn a fusion rule which
minimizes the final loss function. Max-pooling is omitted after
the first convolutional layer in the fusion stage, and a stride
of 2 is used in order to downsample the feature maps while
preserving the spatial information [8]. The use of 3 × 3 filters,
and the absence of max-pooling after the first convolution allows the fusion layer to learn a fusion rule which is somewhat
invariant to spatial mismatches caused by the difference in
imaging modalities. This is due to the fact that the fusion
layer uses 3x3 convolutions to learn relationships between the
features while preserving nearby spatial information. The lack
of max pooling means that these learned spatial relationships
are preserved as not only the maximal response is considered,
while the stride of 2 is used to reduce the feature size. The final
stage of the fusion network consists of two fully connected
layers: the first of which contains 512 channels; while the
second, which performs one-hot binary classification, contains
2 channels.
In a nutshell, the convolutional layers in our network apart
from the fusion layer generally consist of 3 × 3 filters and
follow two rules: 1) The layers with same feature map size
have the same number of filters; and 2) the number of feature
maps increases in the deeper layers, roughly doubling after
each max-pooling layer (except for the last convolutional
stack in each stream). All layers in the network are equipped
with a rectified linear unit (ReLU) as activation function,
except the last fully connected layer, which is activated by
a softmax function. Figure 2 shows the schematic diagram of
the configuration of our network.
B. Loss Function
We
make
use
of
the
binary
cross-entropy
loss
for
training
our
network.
Let
X
=
opt
opt
sar
sar
opt
(xsar
,
x
),
(x
,
x
),
...,
(x
,
x
)
be
a
set
of
SARn
n
1
2
1
2
opt
optical patch pairs, where xsar
∈ RD×D , ∀i = 1, .., n
i , xi

IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, IN PRESS

opt
and yi is the one-hot label for the pair (xsar
i , xi ) (with
[1, 0] denoting a dissimilar pair, and [0, 1] denoting a similar
pair). We then seek to minimize the binary cross-entropy loss
n
1X
opt
(1)
yi · log f (xsar
E=
i , xi , θ)
n i=1
opt
f (xsar
i , xi , θ)

3

SAR Image Stack

(2-D viewing angles)

TomoSAR
Inversion

Stereo
Matching

TomoSAR Point Cloud

Optical Point Cloud

(3-D UTM)

denotes the output vector of the softmax
where
opt
function when comparing the input pair (xsar
i , xi ) with the
current network parameters θ.
C. Configuration Details
Figure 2 depicts the full configuration of our network. Apart
from the previously discussed architecture we also make use of
batch normalization (BN) after the activation function of each
convolutional layer. This leads to an increase in the training
speed and reduces the effects of internal covariate shift. In
order to reduce over-fitting during training, we made use of
L2 -regularization, with λ = 0.001, for the convolution kernels
of the SAR and optical streams, and dropout with a rate of
0.7 for the first fully connected layer.
III. AUTOMATIC PATCH P OOL G ENERATION
For training and testing purposes, a large pool of corresponding and non-corresponding SAR and optical image
patches is needed. While the classical work on deep matching
for optical imagery can usually rely on easy-to-achieve optical
patch pools (see, for example the Phototourism patch dataset
[9], [4]), annotating corresponding patches in VHR optical
and SAR imagery of complex urban scenes is a highly nontrivial task even for experienced human experts. Thus, one of
the contributions of this letter is the introduction of a fully
automatic procedure for SAR-optical patch pool generation.
A. The “SARptical” Framework
In order to solve the challenge of automatic dataset generation, we resort to the so-called “SARptical” framework of
Wang et al. [10], an object-space-based matching procedure
developed for mapping textures from optical images onto 3D
point clouds derived from SAR tomography. The core of
this algorithm is to match the SAR and the optical images
in 3D space in order to deal with the inevitable differences
caused by different geometrical distortions. Usually, this would
require an accurate digital surface model (DSM) of the area
to link homologue image parts via a known object space. In
contrast, the approach in [10] creates two separate 3D point
clouds – one from SAR tomography and one from optical
stereo matching –, which are then registered in 3D space to
form a joint (“SARptical”) point cloud, which serves as the
necessary representation of the object space. The flowchart of
the approach can be seen in Fig. 3. In order to estimate the 3D
positions of the individual pixels in the images, the algorithm
requires an interferometric stack of SAR images, as well as at
least a pair of optical stereo images. The matching of the two
point clouds in 3D guarantees the matching of the SAR and
the optical images. Finally, we can project the SAR image into
the geometry of the optical image via the “SARptical” point
cloud, and vice versa.

Optical Images

(2-D range, azimuth)

(3-D UTM)

3-D
Matching
Matched Point Clouds
(3-D UTM)

Fig. 3. Flowchart of the patch-pool generation procedure.
Fig. 2

B. Data Preparation
For the work presented in this letter, we made use of a stack
of 109 TerraSAR-X high resolution spotlight images of the city
of Berlin, acquired between 2009 and 2013 with about 1 meter
resolution, and of 9 UltraCAM optical images of the same
area with 20cm ground spacing. After the reconstruction of
the “SARptical” 3D point cloud, 8,840 pixels with high SNR
(> 5dB) were uniformly sampled from the non-locally filtered
master SAR amplitude image and projected into the individual
optical images, yielding a total of 10,108 corresponding optical
pixels. The reason for the difference in pixel numbers is that
each of the 9 optical multi-view stereo images is acquired
from a different viewing angle, making it possible for each
SAR image pixel to have up to 9 corresponding optical image
pixels. The actual number of corresponding optical pixels
is dependent on the visibility of the SAR pixel from the
respective optical point of view.
All SAR patches are centered at their corresponding SAR
image pixels. Their size is fixed at 112 × 112 pixels, with a
pixel spacing of about 1 m. In analogy, the optical patches are
centered at the corresponding optical pixels. After resampling
to adjust the pixel spacing, the SAR patches were rotated, so
that both patches align with each other as a first approximation.
In order to reduce bias when training our network we
randomly selected a single correct optical correspondence for
each SAR image patch during the final dataset preparation.
In addition, we randomly assign one wrong optical correspondence to each SAR patch in order to create negative examples.
Thus, eventually, we end up with 17,680 SAR-optical patch
pairs (see Fig. 1 for an example for the class of correct
matches).
As final pre-processing steps, the optical patches were
converted to gray-scale, and all patches were normalized [11]
to a radiometric range of [0; 1] with subsequent subtraction of
their means.
C. Patch Pool Partitioning
In order to provide a fair experimental design, we partition
the patch pool in the following manner: 9724 (55%) of
the patch-pairs are used as training dataset, 2652 (15%) as
validation set, and 5304 (30%) as test dataset. It has to be noted
that we do not partition the patch pool on a purely randomized

IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, IN PRESS

IV. E XPERIMENTS & R ESULTS
A. Training Details
The network was trained using the Adam [12] optimization
algorithm as it is computationally efficient and exhibits faster
convergence than standard stochastic gradient descent methods. The optimization hyper-parameters are fixed to β1 = 0.9,
β2 = 0.999 with a learning rate of αt = 0.0009. The learning
rate was found via a grid search method on our training and
validation data, while the β−parameters were kept at their
recommended values. Prior to training the network weight vectors were initialized using the truncated uniform distribution
described in [13], and the bias vectors were initialized with
zero values. Training was conducted with 2 Nvidia TitanX
GPUs using class balanced, mini-batches of 64 SAR-optical
patch pairs (32 corresponding, 32 non-corresponding pairs)
over 30 epochs; training took on average 25 minutes, with a
single forward pass taking around 3ms to complete.
We trained five versions of our proposed network, each at
a different patch size, in order to evaluate the effect of patch
size on classification accuracy. Patch cropping was done onthe-fly with the new patch being cropped from the center
of a larger patch - this was done as the center pixel is the
point of correspondence between the SAR and optical patch.
Furthermore, we seeded our random number generator with a
fixed value of 0, at the start of training for each patch size, in
order to prevent the randomization effects between networks.

the performance of our proposed network in the presence of
large numbers of potential mismatches.
In Fig. 5a we can see a matrix depicting the similarity scores
of the various pair comparisons, where corresponding SAR
and optical patches are given the same index number. It should
be noted that in determining a binary value for correspondence,
a threshold is applied to these similarity scores. Figure 5b
depicts the sorted scores for all non-similar optical patches,
making it easier to see the number and strength of incorrect
matches in the patch pool.
TABLE I
C ONFUSION MATRIX VALUES FOR DIFFERENT PATCH SIZES

TP
TN

64

76

88

100

112

46.6%
86.2%

61.6%
88.0%

66.0%
88.2%

69.8%
86.0%

82.2%
89.8%

Accuracy Curve for Various Patch Sizes

0.9

64x64
76x76
88x88
100x100
112x112

0.85
0.8
Accuracy

basis but rather resort to a deterministic partitioning method
in order to avoid positively biased test results. The full extent
SAR and optical images are first deterministically partitioned
and then each partition is processed to generate positive and
negative samples for training, validation and testing respectively.

4

0.75
0.7
0.65
0.6
0.55
0.5
0

0.1

0.2

0.3

0.4
0.5
0.6
False Positive Rate

0.7

0.8

0.9

1

Fig. 4. Comparison of different patch sizes.

B. Evaluation Results
We evaluate the proposed network with different input patch
sizes, using our testing patch pool (described in Section III)
which has further been cropped around the center pixel to
produce new testing pools with different patch sizes.
The accuracy vs. false positive rate curves corresponding
to different patch sizes can be seen in Fig. 4. Table I reports
the corresponding confusion matrix values for our proposed
network evaluated with each patch size; it is to be noted that
the confusion matrix is reflective of the network at the point
of highest overall performance for each patch size.

a

b

Fig. 5. Results of key-point matching experiment. a) A confusion matrix
showing the matching scores for all SAR and optical key-point patches. b)
depicts the spread of incorrect matches, ordered by the similarity score.

V. D ISCUSSION
C. Key-Point Matching Results
In order to evaluate the proposed network’s performance
in a real-world, key-point matching scenario we selected 100
neighboring TomoSAR key-points in the SAR image and
extracted the corresponding SAR and optical patch pairs. We
selected these key-points from a localized area within our test
set, so as to reproduce the conditions found in a real-world
key-point matching application. We then compared every SAR
and optical patch in the selected patch set in order to determine

Generally, the results summarized in Section IV-B indicate
a promising discriminative power of the proposed network.
However, the following major points must be considered when
interpreting the results.
A. Influence of the Patch Size
As Tab. I and Fig. 4 clearly indicate, the patch size strongly
affects the discriminative power of the network. This result is
likely due to the effects of distortions in SAR images, which

IEEE GEOSCIENCE AND REMOTE SENSING LETTERS, IN PRESS

True Positives

False Positives

5

False Negatives

True Negatives

Fig. 6. Exemplary patch correspondence results.

are acquired in a range-based imaging geometry. Thus when
patches are cropped to smaller regions we likely crop out
defining features which are used for matching between the
SAR and optical domain. This can be intuitively understood
by referring to Fig. 1 where we can see the effects of layover
and multi-path reflections of the building in the SAR image,
and a near top down view of the same building in the optical
image. Taking away explanatory context will thus render the
matching more difficult. All further discussion will be with
reference to the results we obtained using the largest patch
size, 112 pixels.
B. Comments on the Discriminative Power of the Proposed
Network
In summary, our approach obtains an accuracy exceeding
77% on a separate test dataset when fixing the false positive
rate to 5%, which falls into the same order of magnitude as
what can be achieved using the powerful handcrafted HOPC
descriptor in combination with an L2 -norm cost function [3].
Furthermore, our approach produced a clear diagonal pattern
in Fig. 5a which depicts its ability to accurately determine
the correct correspondence in a key-point matching scenario.
Upon further investigation it was found that the network
achieved 43% top-1 matching accuracy and 74% top-3 accuracy, while 8% of points had no valid matches detected within
the key-point set. This was found to be due to large amounts
of layover and extreme differences in view point between the
SAR and optical patches, see Fig. 6 False Negatives.
C. Possible Reasons for False Predictions
From the randomly chosen prediction examples displayed
in Fig. 6 it can be observed that many of the false positives
and false negatives are erroneously matched under extreme
differences in viewing angle between the SAR and optical
patches. While this may partially be solvable by resorting to
larger patch sizes, providing valuable context, there might be
a need to exclude image parts with all too strong distortions
from further processing.
VI. C ONCLUSION
In this paper, a pseudo-siamese CNN for learning to identify
corresponding patches in very high resolution SAR and optical
images in a fully automatic manner has been presented. A

first evaluation has shown promising potential with respect
to multi-sensor key-point matching procedures. In order to
ensure transferability to other applications not based on keypoints, e.g. dense matching, we will work on the generation of
additional training patches, whose center pixel does not rely on
specific key-points. In addition, we will test the approach on
data coming from a completely different sources. In the end,
we expect our work to help paving the way for generalized
SAR-optical image matching procedures.
ACKNOWLEDGMENT
We thank the NVIDIA Corporation for donating the Titan
X Pascal GPU used in this research.
R EFERENCES
[1] K. Mikolajczyk and C. Schmid, “A performance evaluation of local
descriptors,” IEEE Trans. Pattern Anal. Mach. Intell., vol. 27, no. 10,
pp. 1615–1630, 2005.
[2] M. Schmitt and X. X. Zhu, “Data fusion and remote sensing: An evergrowing relationship,” IEEE Geosci. Remote Sens. Mag., vol. 4, no. 4,
pp. 6–23, 2016.
[3] Y. Ye, J. Shan, L. Bruzzone, and L. Shen, “Robust registration of
multimodal remote sensing images based on structural similarity,” IEEE
Trans. Geosci. Remote Sens., vol. 55, no. 5, pp. 2941–2958, 2017.
[4] X. Han, T. Leung, Y. Jia, R. Sukthankar, and A. C. Berg, “MatchNet:
Unifying feature and metric learning for patch-based matching,” in Proc.
CVPR, Boston, U.S.A., 2015.
[5] S. Zagoruyko and N. Komodakis, “Learning to compare image patches
via convolutional neural networks,” in Proc. CVPR, Boston, U.S.A.,
2015.
[6] L. Mou, M. Schmitt, Y. Wang, and X. Zhu, “A CNN for the identification
of corresponding patches in SAR and optical imagery of urban scenes,”
in Proc. JURSE, Dubai, U.A.E., 2017.
[7] K. Simonyan and A. Zisserman, “Very deep convolutional networks for
large-scale image recognition,” in Proc. ICLR, arXiv:1409.1556, San
Diego, U.S.A., 2014.
[8] A. Handa, M. Bloesch, S. Stent, and A. Davison, “gvnn: Neural network
library for geometric computer vision,” in Proc. ECCV Workshops,
Amsterdam, Netherlands, 2016, pp. 67–82.
[9] S. A. J. Winder and M. Brown, “Learning local image descriptors,” in
Proc. CVPR, Minneapolis, U.S.A., 2007.
[10] Y. Wang, X. Zhu, B. Zeisl, and M. Pollefeys, “Fusing meter-resolution
4-D InSAR point clouds and optical images for semantic urban infrastructure monitoring,” IEEE Trans. Geosci. Remote Sens., vol. 55, no. 1,
pp. 1–13, 2016.
[11] Q. Wang, C. Zou, Y. Yuan, H. Lu, and P. Yan, “Image registration by
normalized mapping,” Neurocomputing, vol. 101, pp. 181–189, 2013.
[12] D. Kingma and J. Ba, “Adam: A method for stochastic optimization,”
in Proc. ICLR, arXiv:1412.6980v9, San Diego, U.S.A., 2015.
[13] X. Glorot and Y. Bengio, “Understanding the difficulty of training deep
feedforward neural networks,” in Proc. AISTATS, Sardinia, Italy, 2010,
pp. 249–256.

