Learning a Distance Metric from Relative Comparisons between Quadruplets of Images

Abstract

This paper is concerned with the problem of learning a distance metric by considering meaningful and discriminative distance constraints in some contexts where rich information between data is provided. Classic metric learning approaches focus on constraints that involve pairs or triplets of images. We propose a general Mahalanobis-like distance metric learning framework that exploits distance constraints over up to four different images. We show how the integration of such constraints can lead to unsupervised or semi-supervised learning tasks in some applications. We also show the benefit on recognition performance of this type of constraints, in rich contexts such as relative attributes, class taxonomies and temporal webpage analysis.

Appendix 1: Solver for the Vector Optimization Problem

We describe here the optimization process when the goal is to learn a dissimilarity function $\mathscr {D}_\mathbf{w }$ parameterized by a vector $\mathbf w$.

1.1 Primal Form of the Optimization Problem

We first rewrite Eq. (14) in the primal form in order to use the efficient and scalable primal Newton method (Chapelle and Keerthi 2010).

The first two constraints of Eq. (14) over $\mathcal {S}$ and $\mathcal {D}$ try to satisfy Eqs. (9) and (10). They are equivalent to $y_{ij}(\mathscr {D}_\mathbf{w } (\mathcal {I}_i, \mathcal {I}_j) - b) \ge 1 - \xi _{ij}$ where $y_{ij} = 1 \Longleftrightarrow (\mathcal {I}_i, \mathcal {I}_j) \in \mathcal {D}$ and $y_{ij} = -1 \Longleftrightarrow (\mathcal {I}_i, \mathcal {I}_j) \in \mathcal {S}$. Equation (14) can then be rewritten equivalently:

$$\begin{aligned} \begin{aligned} \min _{(\mathbf w , b)} \frac{1}{2} (\Vert \mathbf w \Vert _2^2 + b^2)&+ C_p \sum _{(\mathcal {I}_i, \mathcal {I}_j) \in \mathcal {S} \cup \mathcal {D}} L_1(y_{ij},\mathscr {D}_\mathbf{w } (\mathcal {I}_i, \mathcal {I}_j) - b) \\&+ C_q \sum _{q \in \mathcal {N}} L_{\delta _q}(1,\mathscr {D}_\mathbf{w } (\mathcal {I}_k, \mathcal {I}_l)- \mathscr {D}_\mathbf{w } (\mathcal {I}_i, \mathcal {I}_j)) \\ {}&s.t. \quad \mathbf w \in \mathscr {C}^d, b \in \mathscr {C} \end{aligned} \end{aligned}$$

(27)

where $L_1$ and $L_{\delta _q}$ are loss functions and $y_{ij} \in \{ -1 ; 1 \}$. In particular, for Eqs. (14) and (27) to be strictly equivalent, they have to correspond to the classic hinge loss function $L_\delta (y,t) = \max (0,\delta -yt)$. We actually use a differentiable approximation of this function to have good convergence properties (Chapelle 2007; Chapelle and Keerthi 2010).

For convenience, we rewrite some variables:

• $\varvec{\omega } = [\mathbf w ^\top , b]^\top$ is the concatenation of $\mathbf w$ and b in a single $(d+1)$-dimensional vector. We note $e = d+1$ and then have $\varvec{\omega } \in \mathbb {R}^e$.

• $\mathbf c _{ij} = [(\varPsi (\mathcal {I}_{i}, \mathcal {I}_{j}))^\top ,-1]^\top$ is the concatenation vector of $\varPsi (\mathcal {I}_{i}, \mathcal {I}_{j})$ and $-1$. We also have $\mathbf c _{ij} \in \mathbb {R}^e$.

• $p = (\mathcal {I}_{i}, \mathcal {I}_{j}) \Longleftrightarrow \mathbf c _p = \mathbf c _{ij}$ and $y_p = y_{ij}$.

• $q = (\mathcal {I}_i, \mathcal {I}_j, \mathcal {I}_k, \mathcal {I}_l) \Longleftrightarrow \mathbf z _q = \mathbf x _{kl} - \mathbf x _{ij}$.

Equation (27) can be rewritten equivalently with these variables:

$$\begin{aligned} \begin{aligned} \min _{\varvec{\omega } \in \mathscr {C}^e} \frac{1}{2} \Vert \varvec{\omega } \Vert _2^2&+ C_p \sum _{p \in \mathcal {S} \cup \mathcal {D}} L_1(y_{p},\varvec{\omega }^\top \mathbf c _p) \\&+ C_q \sum _{q \in \mathcal {N}} L_{\delta _q}(1, \varvec{\omega }^\top \mathbf z _q) \end{aligned} \end{aligned}$$

(28)

By choosing such a regularization, our scheme may be compared to a RankSVM (Chapelle and Keerthi 2010), with the exception that the loss function $L_{\delta _q}$ works on quadruplets. The complexity of this convex problem w.r.t. $\varvec{\omega }$ is linear in the number of constraints (i.e., the cardinality of $\mathcal {N} \cup \mathcal {D} \cup \mathcal {S}$). It can be solved with a classic or stochastic (sub)gradient descent w.r.t. $\varvec{\omega }$ depending on the number of constraints. The number of parameters to learn is small and grows linearly with the input space dimension, limiting overfitting (Mignon and Jurie 2012). It can also be extended to kernels (Chapelle and Keerthi 2010).

We describe in the following how to apply Newton method (Keerthi and DeCoste 2005; Chapelle 2007; Chapelle and Keerthi 2010) to solve Eq. (28) with good convergence properties. The primal Newton method (Chapelle and Keerthi 2010) is known to be fast for SVM classifier and RankSVM training. As our vector model is an extension of the RankSVM model, the learning is then also fast.

1.2 Loss Functions

Let us first describe loss functions that are appropriate for Newton method. Since the hinge loss function is not differentiable, we use differentiable approximations of $L_1$ and $L_{\delta _q}$ inspired by the Huber loss function.

For simplicity, we also constrain the domain of $\delta _q$ to be 0 or 1 (i.e., $\delta _q \in \{ 0, 1\}$). The set $\mathcal {N}$ can then be partitioned as two sets $\mathcal {A}$ and $\mathcal {B}$ such that for all:

• $q \in \mathcal {N}, \delta _q = 1 \Longleftrightarrow q \in \mathcal {A}$

• $q \in \mathcal {N}, \delta _q = 0 \Longleftrightarrow q \in \mathcal {B}$

In Eq. (28), we consider $t_p = \varvec{\omega }^\top \mathbf c _p$ or $t_q = \varvec{\omega }^\top \mathbf z _q$. Without loss of generality, let us consider $t_r$ with $r \in \varvec{\beta }$ (with $\varvec{\beta } = \mathcal {S}$, $\mathcal {D}$, $\mathcal {A}$ or $\mathcal {B}$) and $y \in \{ -1, +1 \}$. Our loss functions are written:

$$\begin{aligned} L_1^h(y,t_r)= & {} \left\{ \begin{array}{l l l l} 0 &{}\text { if } &{} y t_r > 1 + h &{} \text {set: } \varvec{\beta }_{1,y}^0\\ \frac{(1+h-yt_r)^2}{4h} &{}\text { if } &{} | 1 - y t_r | \le h &{} \text {set: } \varvec{\beta }_{1,y}^Q \\ 1 - y t_r &{} \text { if } &{} y t_r < 1 - h &{} \text {set: } \varvec{\beta }_{1,y}^L\\ \end{array} \right. \end{aligned}$$

(29)

$$\begin{aligned} L_0^h(y,t_r)= & {} \left\{ \begin{array}{l l l l} 0 &{}\text { if } &{} y t_r > 0 &{} \text {set :} \varvec{\beta }_{0,y}^0\\ \frac{t_r^2}{4h} &{}\text { if } &{} | - h - y t_r | \le h &{} \text {set :} \varvec{\beta }_{0,y}^Q\\ - h - y t_r &{} \text { if } &{} y t_r < - 2 h &{} \text {set :} \varvec{\beta }_{0,y}^L\\ \end{array} \right. \nonumber \\ \end{aligned}$$

(30)

where $h \in [0.01,0.5]$. In all our experiments, we set $h = 0.05$.

As described in Chapelle (2007), $L_1^h$ is inspired from the Huber loss function, it is a differentiable approximation of the hinge loss ($L_1(y,t) = \max (0,1-yt)$) when $h \rightarrow 0$. Similarly, $L_0^h$ is a differentiable approximation when $h \rightarrow 0$ of $L_0(y,t) = \max (0,-yt)$, the adaptation of the hinge loss that considers the absence of security margin. Given set $\varvec{\beta }$ and $y \in \{ -1, +1\}$, we can infer three disjoint sets:

• $\varvec{\beta }_{i,y}^0$ is the subset of elements in $\varvec{\beta }$ that have zero loss in $L_i^h(y,\cdot )$.

• $\varvec{\beta }_{i,y}^Q$ is the subset of elements in $\varvec{\beta }$ that are in the quadratic part of $L_i^h(y,\cdot )$.

• $\varvec{\beta }_{i,y}^L$ is the subset of elements in $\varvec{\beta }$ in the non-zero loss linear part of $L_i^h(y,\cdot )$.

1.3 Gradient and Hessian Matrices

By considering $L_1 = L_1^h$ and $L_0 = L_0^h$ in Eq. (28), the gradient $\triangledown \in \mathbb {R}^e$ of Eq. (28) w.r.t. $\varvec{\omega }$ is:

$$\begin{aligned} \begin{aligned} \triangledown&= \varvec{\omega } + \frac{C_p}{2h} \sum _{p \in ({\mathcal {S} \cup \mathcal {D}})_{1,y_p}^Q} (\varvec{\omega }^{\top } \mathbf c _p - (1+h)y_p) \mathbf c _p \\&\quad - C_p \sum _{p \in ({\mathcal {S} \cup \mathcal {D}})_{1,y_p}^L} y_p \mathbf c _p + \frac{C_q}{2h} \sum _{q \in {\mathcal {A}}_{1,1}^Q} (\varvec{\omega }^{\top } \mathbf z _q - (1+h)) \mathbf z _q \\&\quad + \frac{C_q}{2h} \sum _{q \in {\mathcal {B}}_{0,1}^Q} (\varvec{\omega }^{\top } \mathbf z _q) \mathbf z _q - C_q \sum _{q \in (\mathcal {A}_{1,1}^L \cup \mathcal {B}_{0,1}^L)} \mathbf z _q \end{aligned} \end{aligned}$$

(31)

and the Hessian matrix $\mathbf H \in \mathbb {R}^{e \times e}$ of Eq. 28 w.r.t. $\varvec{\omega }$ is:

$$\begin{aligned} \mathbf H = \mathbf I _e + \frac{C_p}{2h} \sum _{p \in ({\mathcal {S} \cup \mathcal {D}})_{1,y_p}^Q} \mathbf c _p \mathbf c _p^{\top } + \frac{C_q}{2h} \sum _{q \in (\mathcal {A}_{1,1}^Q \cup \mathcal {B}_{0,1}^Q)} \mathbf z _q \mathbf z _q^{\top } \end{aligned}$$

(32)

where $\mathbf I _e \in \mathbb {R}^{e \times e}$ is the identity matrix. $\mathbf H$ is the sum of a positive definite matrix ($\mathbf I _e$) and of positive semi-definite matrices. $\mathbf H$ is then positive definite, and thus invertible (because every positive definite matrix is invertible).

Proof

$\mathbf H$ can be written $\mathbf H = \mathbf I _e + \mathbf B$ with $\mathbf B \in \mathbb {R}^{e \times e}$ a positive semi-definite matrix. For all vector $\mathbf z \in \mathbb {R}^{e}$, we have $\mathbf z ^{\top } \mathbf H \mathbf z = \mathbf z ^{\top }{} \mathbf I _e\mathbf z + \mathbf z ^{\top } \mathbf B \mathbf z$. By definition of positive (semi-)definiteness, we have the following property: for all nonzero $\mathbf z \in \mathbb {R}^{e}$, $\mathbf z ^{\top }{} \mathbf I _e\mathbf z > 0$ and $\mathbf z ^{\top } \mathbf B \mathbf z \ge 0$. Then for all nonzero $\mathbf z \in \mathbb {R}^{e}$, $\mathbf z ^{\top } \mathbf H \mathbf z > 0$. $\mathbf H$ is then a positive definite matrix. $\square$

The global learning scheme is described in Algorithm 2. The step size $\eta _t > 0$ can be set to 1 and unchanged as in Chapelle (2007), or optimized at each iteration through line search (see Sect. 9.5.2 in Boyd and Vandenberghe 2004). The parameter $\epsilon \ge 0$ determines the stopping criterion by controlling the $\ell _2$-norm of the difference of $\varvec{\omega }$ between iteration t and $t-1$.

Complexity Computing the Hessian takes $O(\sigma e^2)$ time (where $\sigma = |({\mathcal {S} \cup \mathcal {D}})_{1,y_p}^Q|+|(\mathcal {A}_{1,1}^Q \cup \mathcal {B}_{0,1}^Q)|$) and solving the linear system is $O(e^3)$ because of the inversion of $\mathbf H _t \in \mathbb {R}^{e \times e}$. This can be prohibitive if e is large but we restrict $e \le 1001$ in our experiments; the inversion of $\mathbf H _t$ is then very fast. Other optimization methods are proposed in Chapelle and Keerthi (2010) (e.g., a truncated Newton method) if e is large.

It can be noticed that Newton method is appropriate for unconstrained problems, where the inclusion of $\mathbf H ^{-1}$ at each iteration allows to converge faster to the global minimum. When $\mathscr {C}^e$ is $\mathbb {R}_+^e$, Eq. (28) is a constrained problem and the minimum of the unconstrained problem is not necessarily the minimum of the constrained problem. In Eq. (28), since our loss functions are linear almost everywhere on their domain, the Hessian of the problem is close to the identity matrix and it is affected almost exclusively by the regularization term. This is why applying a projected Newton method is not a major issue in our case. If computing the inverse of the Hessian is too much expensive, the Hessian can be omitted and a classic projected gradient method can be used.