Integral Geometric Dual Distributions of Multilinear Models

TL;DR

This paper introduces an integral geometric method to compute dual distributions of parameter estimates in multilinear models, enabling analysis of uncertainty and confidence intervals in feature distributions.

Contribution

It presents a novel integral geometric approach using Radon transforms to derive analytical dual distributions for various multilinear models.

Findings

Allows computation of dual distributions from parameter estimates

Enables analysis of uncertainty in feature and training data distributions

Provides a framework for deriving confidence intervals of model parameters

Abstract

We propose an integral geometric approach for computing dual distributions for the parameter distributions of multilinear models. The dual distributions can be computed from, for example, the parameter distributions of conics, multiple view tensors, homographies, or as simple entities as points, lines, and planes. The dual distributions have analytical forms that follow from the asymptotic normality property of the maximum likelihood estimator and an application of integral transforms, fundamentally the generalised Radon transforms, on the probability density of the parameters. The approach allows us, for instance, to look at the uncertainty distributions in feature distributions, which are essentially tied to the distribution of training data, and helps us to derive conditional distributions for interesting variables and characterise confidence intervals of the estimates.

Figures

Tables

Table 1. TABLE I: The trilinearities, adopted from [ 1 ] .

Correspondence	Relation		dof
three points	$x^{i}x_{}^{\prime }{}_{}{}^{j}x_{}^{\prime \prime }{}_{}{}^k{ϵ}_{jqs}{ϵ}_{krt}{𝒯}_i^{qr}$=	$0_{st}$	4
two points, one line	$x^{i}x_{}^{\prime }{}_{}{}^j{l}_r^{\prime \prime }{ϵ}_{jqs}{𝒯}_i^{qr}$=	$0_s$	2
one point, two lines	$x^i{l}_q^{\prime }{l}_r^{\prime \prime }{𝒯}_i^{qr}$=	$0$	1
three lines	$l_p{l}_q^{\prime }{l}_r^{\prime \prime }{ϵ}^{piw}{𝒯}_i^{qr}$=	$0^w$	2

Equations

$$f(\dots ,\alpha \mathbf{x}+\beta \mathbf{y},\dots )=\alpha f(\dots ,\mathbf{x},\dots )+\beta f(\dots ,\mathbf{y}\dots ),$$

$$f({\mathbf{x}}_1,{\mathbf{x}}_2,\dots ,{\mathbf{x}}_n;\theta )=0,$$

$$\mathbf{y}\widehat{=}{\mathbf{x}}_1\otimes {\mathbf{x}}_2\otimes \dots \otimes {\mathbf{x}}_n$$

$${\theta }^{\mathrm{T}}{\mathbf{T}}_l\mathbf{y}\equiv {\theta }^{\mathrm{T}}{\mathbf{y}}_l=0,$$

$$p(\theta )\propto \exp \left(-\frac{1}{2}{\theta }^{\mathrm{T}}{\mathbf{C}}_{\theta }^{†}\theta \right),$$

$${\theta }^{\mathrm{T}}\mathbf{y}=0,$$

$$\begin{array}{rl}{\theta }^{\mathrm{T}}{\mathbf{C}}_{\theta }^{†}\theta & ={\theta }^{\mathrm{T}}\mathbf{U}{\mathbf{\Lambda }}^{†}{\mathbf{U}}^{\mathrm{T}}\theta \\ & =\underset{\triangleq {\theta }^{\prime \mathrm{T}}}{\underbrace{{\theta }^{\mathrm{T}}\tilde{\mathbf{U}}\left(\begin{array}{cc}{\tilde{\mathbf{\Lambda }}}^{-1/2}& 0\\ 0^{\mathrm{T}}& 1\end{array}\right)}}\underset{\triangleq \tilde{\mathbf{I}}}{\underbrace{\left(\begin{array}{cc}\mathbf{I}& 0\\ 0^{\mathrm{T}}& 0\end{array}\right)}}\underset{\triangleq {\theta }^{\prime }}{\underbrace{\left(\begin{array}{cc}{\tilde{\mathbf{\Lambda }}}^{-1/2}& 0\\ 0^{\mathrm{T}}& 1\end{array}\right){\tilde{\mathbf{U}}}^{\mathrm{T}}\theta }}\\ & ={\theta }^{\prime \mathrm{T}}\tilde{\mathbf{I}}{\theta }^{\prime },\end{array}$$

$$\begin{array}{rl}0& ={\theta }^{\mathrm{T}}\mathbf{y}=\underset{{\theta }^{\prime \mathrm{T}}}{\underbrace{{\theta }^{\mathrm{T}}\tilde{\mathbf{U}}\left(\begin{array}{cc}{\tilde{\mathbf{\Lambda }}}^{-1/2}& 0\\ 0^{\mathrm{T}}& 1\end{array}\right)}}\underset{\triangleq {\mathbf{y}}^{\prime }}{\underbrace{\left(\begin{array}{cc}{\tilde{\mathbf{\Lambda }}}^{1/2}& 0\\ 0^{\mathrm{T}}& 1\end{array}\right){\tilde{\mathbf{U}}}^{\mathrm{T}}\mathbf{y}}}={\theta }^{\prime \mathrm{T}}{\mathbf{y}}^{\prime },\end{array}$$

$$\begin{array}{rl}p(\rho ,\varphi )& =\frac{\Gamma (M/2-1/2){\mathrm{e}}^{-\frac{1}{2}{\rho }^{-2}}{\prod }_{i=1}^{M-3}{\sin }^{M-2-i}({\varphi }_i)}{\sqrt{2{\pi }^M}{\rho }^2},\end{array}$$

$$s=-\frac{y_N^{\prime }}{\mathrm{sign}(y_{N-1}^{\prime })\sqrt{{\mathbf{y}}^{\prime \mathrm{T}}{\mathbf{y}}^{\prime }-{y^{\prime }}_N^2}}.$$

$$s=-\frac{1}{\rho }.$$

$$p(s|\varphi )={\int }_{{{\mathbf{y}}_s^{\prime }}^{\mathrm{T}}{\theta }^{\prime }=0}{p}_{\mathrm{G}}({\theta }^{\prime };{\mathbf{e}}_M,\tilde{\mathbf{I}})\mathrm{d}S=p(s)$$

$$\begin{array}{rl}p(\rho |\varphi )& =\left|\frac{\partial s}{\partial \rho }\right|p\left(s(\rho )|\varphi \right)=\frac{1}{\sqrt{2\pi }{\rho }^2}\exp \left(-\frac{1}{2}{\rho }^{-2}\right).\end{array}$$

$$\begin{array}{r}p(\varphi )=\frac{\Gamma \left(M/2-1/2\right)}{{\pi }^{(M-1)/2}}\prod _{i=1}^{M-3}{\sin }^{M-2-i}({\varphi }_i).\end{array}$$

$$\begin{array}{rl}p(\rho ,\varphi )& =p(\rho |\varphi )p(\varphi ),\end{array}$$

$${\theta }^{\mathrm{T}}{\mathbf{y}}_l=0,l=1,\dots ,L.$$

$$\begin{array}{rl}0& ={\theta }^{\mathrm{T}}\mathbf{Y}=\underset{{\theta }^{\prime \mathrm{T}}}{\underbrace{{\theta }^{\mathrm{T}}\tilde{\mathbf{U}}\left(\begin{array}{cc}{\tilde{\mathbf{\Lambda }}}^{-1/2}& 0\\ 0^{\mathrm{T}}& 1\end{array}\right)}}\underset{\triangleq {\mathbf{Y}}^{\prime }}{\underbrace{\left(\begin{array}{cc}{\tilde{\mathbf{\Lambda }}}^{1/2}& 0\\ 0^{\mathrm{T}}& 1\end{array}\right){\tilde{\mathbf{U}}}^{\mathrm{T}}\mathbf{Y}}}\\ & ={\theta }^{\prime \mathrm{T}}{\mathbf{Y}}^{\prime }.\end{array}$$

$${\theta }_{\min }^{\prime }\triangleq \left(\begin{array}{c}\mathbf{w}\\ 1\end{array}\right)=\left(\begin{array}{c}{\mathbf{A}}^{†}\mathbf{b}\\ 1\end{array}\right).$$

$$\mathbf{P}=\mathbf{I}-{\mathbf{A}}^{†}\mathbf{A},$$

$$\mathbf{P}(\Phi)=\left\{\begin{tabular}[]{ll}$\mathbf{L}(\Phi)\mathbf{L}(\Phi)^{\mathrm{T}}$,&if $K\leq M/2$\\ $\mathbf{I}-\mathbf{L}(\Phi)\mathbf{L}(\Phi)^{\mathrm{T}}$,&otherwise.\end{tabular}\right.$$

$$\begin{array}{rl}p(\mathbf{s},\Phi )& =\frac{{\mathrm{e}}^{-\frac{1}{2}{\mathbf{s}}^{\mathrm{T}}\mathbf{s}}}{2^{(M-K-1)/2}{\pi }^{(M-K+\tilde{K}M-{\tilde{K}}^2-1)/2}}\\ & \times \prod _{k=1}^{\tilde{K}}\Gamma \left((M-2k+1)/2\right)\prod _{i=1}^{M-2k-1}{\sin }^{M-2k-i}({\varphi }_i^k).\end{array}$$

$$\begin{array}{rl}p(\mathbf{s}|\Phi )=& {\int }_V{p}_{\mathrm{G}}({\theta }^{\prime };{\mathbf{e}}_M,\tilde{\mathbf{I}})\mathrm{d}S\\ \text{mbox}=& p(\mathbf{s}),\end{array}$$

$$\begin{array}{rl}p(\Phi )& =p({\Phi }_1,{\Phi }_2,\dots ,{\Phi }_{\tilde{K}})\\ & =p({\Phi }_1)p({\Phi }_2|{\Phi }_1)p({\Phi }_3|{\Phi }_1,{\Phi }_2)\cdots \\ & p({\Phi }_{\tilde{K}}|{\Phi }_1,{\Phi }_2,\dots ,{\Phi }_{\tilde{K}-1}),\end{array}$$

$$\begin{array}{rl}p({\Phi }_k|& {\Phi }_1,{\Phi }_2,\dots ,{\Phi }_{k-1})\\ & =\frac{\Gamma \left(\left(M-2k+1\right)/2\right)}{{\pi }^{(M-2k+1)/2}}\prod _{i=1}^{M-2k-1}{\sin }^{M-2k-i}({\varphi }_i^k),\end{array}$$

$$\begin{array}{rl}p(W)& =p(\mathbf{s},\Phi )=p(\mathbf{s}|\Phi )p(\Phi ),\end{array}$$

$$\begin{array}{rl}p(\mathbf{w}|{\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_K)=& {\int }_W{p}_{\mathrm{G}}({\theta }^{\prime };{\mathbf{e}}_M,\tilde{\mathbf{I}})\mathrm{d}S\\ \text{mbox}=& p(\mathbf{w}),\end{array}$$

$$\mathbf{s}={\tilde{\mathbf{V}}}^{\mathrm{T}}\mathbf{w}\iff \mathbf{w}={\tilde{\mathbf{V}}}^{-\mathrm{T}}\mathbf{s},$$

$$\begin{array}{rl}p(\mathbf{s}|{\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_K)& =|\det \mathrm{J}\mathbf{w}(\mathbf{s})|p(\mathbf{w}(\mathbf{s})|{\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_K)\\ & =\frac{{\mathrm{e}}^{-\frac{1}{2}\Vert {\tilde{\mathbf{V}}}^{-\mathrm{T}}\mathbf{s}{\Vert }^2}}{(2\pi {)}^{K/2}|\det (\tilde{\mathbf{V}})|}.\end{array}$$

$$\begin{array}{rl}p({\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_K)& =\frac{\Gamma (M/2-1/2{)}^K}{{\pi }^{K(M-1)/2}}.\end{array}$$

$$\begin{array}{rl}p(\mathbf{s},{\mathbf{v}}_1,& {\mathbf{v}}_2,\dots ,{\mathbf{v}}_K)\\ =& p(\mathbf{s}|{\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_K)p({\mathbf{v}}_1,{\mathbf{v}}_2,\dots ,{\mathbf{v}}_K)\\ =& \frac{\Gamma (M/2-1/2{)}^K{\mathrm{e}}^{-\frac{1}{2}\Vert {\tilde{\mathbf{V}}}^{-\mathrm{T}}\mathbf{s}{\Vert }^2}}{\sqrt{2^K{\pi }^{KM}}|\det (\tilde{\mathbf{V}})|}.\end{array}$$

Full text

Keywords: Duality, Radon Transform, Uncertainty, Confidence Intervals, Integral Geometry, Computer Vision.

1 Introduction

An essential part of geometric computer vision are multilinear models which include e.g. homographies, multiple view tensors, quadric surfaces, as well as the simple entities of points, lines and planes. There are numerous works related to the estimation of these kinds of multilinear relations, see e.g. [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. When a statistical approach is selected for the estimation, one should be able to compute the geometric model parameters along with their uncertainty distribution. This paper considers how these parameter distributions of multilinear geometric entities can be dualised. The simplest form of this dualisation is the transformation of the line-probability-density into a point-probability-density as proposed in [12]. This paper generalises the dualisation approach for general multilinear models. An early, conference version of this paper is [11].

By the way of an example, consider fitting a conic section to a set of points using maximum likelihood estimation. We would be interested in the confidence intervals of the conic, but the normal distribution assumption can be made only for the MLE in the parameter space. However, we will show that by the dualisation of the parameter distribution we will obtain a distribution of points that can be used to plot the selected confidence intervals for the estimated conic section. As the second example, let us consider the trifocal point transfer. Given a point match in two views, it would be interesting to compute the conditional position distribution of the transferred point in the third view if we had the uncertainty information of the trifocal tensor available. In fact, by dualising the trifocal tensor parameter distribution, an exact form for this conditional distribution can be computed, as will be shown later in this paper.

Our approach is closely related to the branch of integral geometry in mathematics. There are two main schools of integral geometry of which the traditional is that of Santaló and Blaschke [13]. The classical example is that the length of a plane curve is the probability of random lines intersecting it. The more recent meaning is the school of Gelfand [14]. It studies integral transforms, modelled with the Radon transform, that relates the underlying geometrical incidence relations by incidence graphs. Our approach seems to be somewhat in between these two schools as we compute (generalised) Radon transforms for the probability densities in such a way that the probability measure is preserved. The dualisation is constructed from the fact that the probability of an element (e.g. a point) is the total probability of all the geometric entities (e.g. planes) that coincide with the element. Then by integrating the distribution of the entity over the affine subspaces corresponding to the selected incidence relation, Radon like integral transforms follow.

As the principal assumption we use the normal distribution assumption for the multilinear model parameters. This is reasonable due to the asymptotic normality property of the maximum likelihood estimator, i.e., due to the fact that, with certain general regularity conditions, the distribution of the maximum likelihood parameter estimator converges in distribution to the normal distribution with the (pseudo)inverse of the Fisher information as the parameter covariance matrix [15, 2]. This makes the approach taken here fundamentally different from the work in [16] where an algebraic linear system and Gaussian approximation in the feature space were used. In contrast to the work in [16], the dual distributions considered here have analytic forms and are exact with the assumptions above.

This paper is organised as follows. In Section 2, we introduce the multilinear models considered in this paper. In Section 4, we derive the dual distributions by first assuming a single constraint equation and then generalise the approach for multiple constraint equations. In Section 5, we show how interesting conditional distributions can be extracted. In Section 6, we compute confidence intervals for conics and compute the point transfer density from two views into the third view. Conclusions are in Section 8.

2 Multilinear Model

We first need to define what we mean by a multilinear model.

Definition 2.1

A function $f:V_{1}\times\cdots\times V_{k}\rightarrow W$, where $V_{1}\ldots,V_{k}$ and $W$ are real vector spaces, is k-linear if it is linear in each of its $k$ arguments:

[TABLE]

for all $\alpha,\beta\in\mathbb{R}$ and $\mathbf{x},\mathbf{y}\in V$.

Definition 2.2

Let $f$ be a $(n+1)$-linear function $f:\mathbb{P}^{m}\times\ldots\times\mathbb{P}^{m}\times\mathbb{P}^{N-1}\rightarrow\mathbb{R}^{L}$. The multilinear model is defined as the relation

[TABLE]

*where $\mathbf{x}_{i}$ are feature vectors and $\theta$ is the parameter vector, where $\mathbf{x}_{i}\in\mathbb{P}^{m}$, $\theta\in\mathbb{P}^{N-1}$ and $i=1,2,\ldots,n$. *

The definition is to be taken in the general sense so that any of the $n$ first arguments may be repeated arbitrary many times so that, e.g., quadratic forms are included.

With a fixed $\theta$, the multilinear model defines a multilinear equation system, where each equation is equivalent to a linear subspace, with co-dimension one, in the space of the joint feature vector.

Definition 2.3

The joint feature vector is defined as the vector $\mathbf{y}\in\mathbb{P}^{N-1}$, such that

[TABLE]

*containing the elements of the tensor product, up to scale, where the repeating elements have been dropped; $\hat{=}$ denotes the correspondence between the expressions, $\otimes$ is the tensor product, and $\mathbf{x}_{i}\in\mathbb{P}^{m}$, $i=1,2,\ldots,n$. *

For the case where the feature vectors are distinct, the mapping from the feature vectors to the joint feature vector is known as Segre embedding.

Without a loss of generality, we assume that each equation $l$ of the multilinear system is in the form of

[TABLE]

where $\mathbf{T}_{l}$ is a matrix defined by the multilinear relation. For instance, the well known point and line incidence relations [1, 17] characterising multiple projective views of a scene can be written in this form.

3 Statistical Model

Let us assume that we have the maximum likelihood estimate $\theta_{0}$, or the corresponding robust estimate [4], for the parameter vector and its the covariance matrix $\mathbf{C}_{\theta}$ available, where $\theta_{0}$ is constrained to lie on the unit hypersphere $\mathbb{S}^{N-1}$ in $\mathbb{R}^{N}$. Using the asymptotic normality property of the MLE, we assume that $\theta\sim N(\theta_{0},\mathbf{C}_{\theta})$. Since $\theta_{0}\in\mathrm{Ker}\ \{\mathbf{C}_{\theta}\}$, $\theta$ has the density function

[TABLE]

where $\cdot^{\dagger}$ is the Moore–Penrose pseudoinverse. That is, we use the tangential, normal approximation for the variation of $\theta$ on the unit hypersphere at $\theta_{0}$.

4 Duals of the Parameter Distributions

In this section, we dualise the normal parameter uncertainty distributions starting from the models defined by a single constraint equation and generalise the point–line case [12] to the point–hyperplane duality (Section 4.1). The general case of multiple constraint equations is considered in Section 4.2.

4.1 Single Constraint Equation

Let $\theta$ be a parameter vector that defines the multilinear model by a single constraint equation

[TABLE]

with the joint feature vector $\mathbf{y}\in\mathbb{P}^{N-1}$.

Our intention is to compute the dual of $p(\theta)$. The duality between the parameter vector $\theta$ and the feature vector $\mathbf{y}$ is illustrated in Fig. 1. As (6) represents the hyperplane $\theta$ in the space of $\mathbf{y}$, it can also be seen as the hyperplane $\mathbf{y}$ in the dual space or space of $\theta$. Now, as the dual distribution of $\theta$, which is a distribution for $\mathbf{y}$, we will identify the total probability (density) of the all the models $\theta$ lying on the hyperplane $\mathbf{y}$ in the dual space. To construct the dual pdf $p(\mathbf{y})$, we first make the whitening transformation for variables. Then we compute the Radon transform, i.e., transform the whitened dual domain by computing the integrals of $p(\theta)$ over all the hyperplanes in the dual space. The subsets of these integrals corresponding to parallel hyperplanes form the conditional probability density (marginal density) conditioned on a fixed normal direction of the hyperplane. Multiplying this conditional density with the pdf of normal directions, which is the uniform density on the half of the unit hypersphere with the whitened Gaussian model, we obtain a valid probability density which can be interpreted in the parameter space of $\mathbf{y}$. The details are below.

{degenerate1}

Let us make the whitening transform from the eigenvalue decomposition of $\mathbf{C}_{\theta}$, which has been constructed so that $\theta_{0}$ corresponds to the last eigenvector. We may write

[TABLE]

where the diagonal matrix $\tilde{\mathbf{\Lambda}}$ contains the $M-1$ non-zero eigenvalues of $\mathbf{C}_{\theta}$, sorted in the descending order, and $\tilde{\mathbf{U}}$ contains the corresponding eigenvectors and the eigenvector representing $\theta_{0}$. Now, $\theta^{\prime}\sim N(\mathbf{e}_{M},\tilde{\mathbf{I}})$, where $\mathbf{e}_{M}$ is the standard basis vector $(0,\ldots,0,1)\in\mathbb{R}^{M}$. Furthermore, for all the feature vectors that are consistent with the model $\theta$, lying on the tangent space, we may write

[TABLE]

where $\mathbf{y}^{\prime}\in\mathbb{P}^{M-1}$ represents a reduced joint feature vector. For simplicity, we now investigate the reduced, transformed model $\theta^{\prime}$. Now we are ready to state our main theorem for the case of the single constraint equation.

Theorem 4.1

Let $\mathbf{y},\theta\in\mathbb{P}^{N-1}$, normalised so that the joint feature vector $\mathbf{y}$ has homogeneous scaling of unity and $\theta$ lies in the tangent space $T_{\theta_{0}}(\mathbb{S}^{N-1})$. Moreover, let $\theta\sim N(\theta_{0},\mathbf{C})$ with the probability density function $p(\theta)$, where the eigenvectors of $\mathbf{C}$ corresponding to the $M-1$ non-zero eigenvalues span $T_{\theta_{0}}$, and $\mathbf{y}^{\prime}$ is the reduced joint feature vector corresponding to $\mathbf{y}$. The dual distribution of $p(\theta)$ has the analytic form

[TABLE]

where the reduced joint feature vector $\mathbf{y}^{\prime}$ is parameterized by the modified hyperspherical coordinates $(\rho,\phi)$ in $\mathbb{R}^{M-1}$.

Proof:

On the basis of the construction above, the variation of $\theta^{\prime}$ occurs only in the $M-1$ dimensional affine subspace $\pi^{\prime}$ perpendicular to $\mathbf{e}_{M}$ in the dual space. Since every point on this tangent hyperplane can be identified with a unique one-dimensional linear subspace of $\mathbb{R}^{M}$, we may regard the tangent hyperplane as a projective space $\mathbb{P}^{M-1}$. Moreover, the points on the tangent plane can be considered to be already in the homogeneous form. Hence, in the following we assume that $\theta^{\prime}\in\mathbb{P}^{M-1}$.

In the dual space, the $M-2$ dimensional hyperplanes $\mathbf{y}^{\prime\mathrm{T}}\theta^{\prime}=0$ embedded in the $M-1$ dimensional affine subspace $\pi^{\prime}\hat{=}\mathbb{R}^{M-1}$, may be parameterised by the signed distance $s$ from the origin and by the unit vector $\mathbf{v}\in\mathbb{S}^{M-2}$, as Fig 2(b) illustrates. The origin can be represented by $\mathbf{e}_{M}\in\mathbb{P}^{M-1}$ and we assume that $\mathbf{v}=\mathbf{v}(\phi)$ where $\phi$ parameterises the normal direction of the $M-2$ dimensional hyperplane $\mathbf{y}^{\prime}$ in $\mathbb{R}^{M-1}$. We choose the sign of $s$ to be equal to sign of the intercept of the hyperplane in the dual space, hence,

[TABLE]

On the other hand, in the reduced joint feature space, the homogeneous vector $\mathbf{y}^{\prime}$ can be parameterised by the signed distance $\rho$ from the origin $\mathbf{e}_{M}$ and the direction $\mathbf{v}$ (Fig 2(a)). The sign of $\rho$ is identified as the last variable sign of the inhomogeneous representation of $\mathbf{y}^{\prime}$. Then we have

[TABLE]

As we assume that a $\theta^{\prime}\sim N(\mathbf{e}_{M},\tilde{\mathbf{I}})$, the probability of the hyper plane ${\mathbf{y}^{\prime}_{s}}^{\mathrm{T}}\theta^{\prime}=0$, conditioned on the direction $\mathbf{v}(\phi)$ in the dual space, is simply the marginal probability of the Gaussian over the hyperplane, i.e.,

[TABLE]

which is a mean zero, 1-D Gaussian with unity variance [12], where $\mathrm{d}S$ denotes the volume differential in the hyperplane.

By using the modified spherical coordinates $\mathbf{y}^{\prime}=\mathbf{y}^{\prime}(\rho,\phi)$ in the reduced joint feature space {automatic} (see Appendix A. Modified Spherical Coordinates in $\mathbb{R}^{N}$)

{manual} (see Appendix A)

and the fact that $s=s(\rho)$, we get

[TABLE]

On the other hand, since we have an isotropic Gaussian distribution, the distribution of normal directions is uniform over a half of the unit sphere. By using the modified spherical coordinates, and marginalising the $M-1$-dimensional Gaussian over the signed radial parameter, we obtain

[TABLE]

Hence,

[TABLE]

and the claim follows. ∎

We have thus derived an analytic form for the probability density of the reduced joint feature vector $\mathbf{y}^{\prime}=\mathbf{y}^{\prime}(\rho,\phi)$ assuming the model (6) and the Gaussian distribution for the parameter vector $\theta$.

4.2 Multiple Constraint Equations

Now we generalise the computations above for multiple constraint equations. Let $\theta$ be a parameter vector of a multilinear model. We define that the set vectors $\mathbf{y}_{l},\ l=1,\ldots,L$ are consistent with the model if and only if

[TABLE]

Hence, $(\mathbf{y}_{1}\ \dots\ \mathbf{y}_{L})^{\mathrm{T}}\theta=\mathbf{Y}^{\mathrm{T}}\theta=\mathbf{0}$, i.e., $\theta\in\mathrm{Ker}\ \{\mathbf{Y}^{\mathrm{T}}\}$ (c.f. the nullspace representation of lines in $\mathbb{R}^{3}$). Similarly as in (8), we define the transformed model $\theta^{\prime}$ and the reduced coefficient matrix $\mathbf{Y}^{\prime}$ so that

[TABLE]

Let $W$ be the affine subspace corresponding to the set $\{\theta^{\prime}\in\mathbb{P}^{M-1}|\mathbf{Y}^{\prime\mathrm{T}}\theta^{\prime}=0\}$ or $W=\{\tilde{\theta}^{\prime}\in\mathbb{R}^{M-1}|\mathbf{A}\tilde{\theta}^{\prime}=\mathbf{b}\}$, where $\theta^{\prime}\hat{=}(\tilde{\theta}^{\prime},1)$ and ${\mathbf{Y}^{\prime}}^{\mathrm{T}}=(\mathbf{A}\ \ -\mathbf{b})$. Assuming that the $W$ contains finite points, it is easy to show that the nearest point to the origin $\mathbf{e}_{M}$ is

[TABLE]

In other words, the affine subspace $W=\mathbf{w}+V$ where $V$ denotes the linear subspace parallel to $W$ and $\dim(V)\equiv K$.

We are interested in the probability density of $\mathrm{Ker}\ \{\mathbf{Y}^{\prime\mathrm{T}}\}$, that is, the total probability (density) of the affine subspace $W$ as it represents the set of all the models $\theta^{\prime}$ that are consistent with $\mathbf{Y}^{\prime}$. To construct this dual density, instead of integrating over all hyperplanes as in the previous subsection, we must generalise the Radon transform by integrating $p(\theta^{\prime})$ over all the $K$-dimensional affine subspaces. To perform this kind of integral transform and to interpret its result as probability densities, we must create a unique parameterisation for affine subspaces. A unique parameterisation could be constructed by the Grassmannian coordinates [18, 16]—it could give certain algebraic benefits that are to be investigated in future.

In this paper, we study Gaussians in the reduced coordinate frame that suggests a simple parameterisation, generalising the result on the single coordinate equation. We represent an affine subspace by the offset vector $\mathbf{w}$ and parameters of the parallel linear subspace $V$. To parameterise $V$, we use the fact that the related orthogonal projection matrix

[TABLE]

projecting onto the subspace is unique. It is well known that a matrix represents an orthogonal projection onto a linear subspace if and only if it is idempotent, $\mathbf{P}^{2}=\mathbf{P}$, and self-adjoint, $\mathbf{P}^{\mathrm{T}}=\mathbf{P}$. Now, to uniquely parameterise the matrix $\mathbf{P}$ we need the following lemma.

Lemma 4.2

*Let $\mathbf{P}$ be an idempotent and symmetric $(M-1)\times(M-1)$ matrix. If $\mathbf{p}_{i}=\mathbf{P}\mathbf{e}_{i}$, $i=1,\ldots,K$ are linearly independent and span the range of $\mathbf{P}$, where $K\leq\frac{M}{2}$, is the dimension of the range, then there is a unique lower triangular $K\times(M-1)$ matrix $\mathbf{L}$ with orthonormal columns and strictly positive diagonal so that $\mathbf{P}=\mathbf{L}\mathbf{L}^{\mathrm{T}}$. *

The proof is in {automatic} Appendix B. Proof of Lemma 1.

{manual} Appendix B.

The lemma suggests that we may parameterise the orthogonal projection matrix $\mathbf{P}$ by parameterising the elements of $\mathbf{L}=\mathbf{L}(\Phi)$ when we may define

[TABLE]

The matrix $\mathbf{L}$ can be parameterised by the modified spherical coordinates {automatic} (see Appendix A. Modified Spherical Coordinates in $\mathbb{R}^{N}$)

{manual} (see Appendix A)

with the radial parameter equal to unity. The first column vector is on the half of the unit sphere $\mathbb{S}^{M-2}$, the second column vector is orthogonal to the first and it has one element less and is hence on $\mathbb{S}^{M-4}$. To parameterise $\mathbb{S}^{M-4}$, we form an orthogonal basis in the orthogonal complement of the first column vector, by creating the orthogonal projection matrix that projects onto the orthogonal complement, and employ the Gram–Schmidt orthonormalisation procedure, as shown in Appendix B, and form the unit sphere in the subspace spanned by the orthonormal vectors. The other columns $l\leq\tilde{K}$ can be similarly parameterised on the half of the unit sphere $\mathbb{S}^{M-2l}$, whereas $\tilde{K}=K$, if $K\leq M/2$, and $\tilde{K}=M-K-1$, otherwise. Let $\Phi_{k}$ denote the vector of modified spherical coordinates of the column $k$ in $\mathbf{L}$. We collect the parameters in the vector $\Phi=(\Phi_{1},\Phi_{2},\ldots,\Phi_{\tilde{K}})$. In total, there are $K^{\prime}=(M-1)\tilde{K}-\tilde{K}^{2}$ free parameters in $\mathbf{L}$.

Now, we are ready to state our theorem about the probability density for the affine subspace $W$.

Theorem 4.3

Let $\mathbf{Y}\in\mathbb{R}^{N\times L}$ and $\theta\in\mathbb{P}^{N-1}$, normalised so that $\theta$ lies in the tangent space $T_{\theta_{0}}(\mathbb{S}^{N-1})$. Moreover, let $\theta\sim N(\theta_{0},\mathbf{C})$ with the probability density function $p(\theta)$, where the eigenvectors of $\mathbf{C}$ corresponding to the $M-1$ non-zero eigenvalues span $T_{\theta_{0}}$, and $\mathbf{Y}^{\prime}$ is the reduced matrix corresponding to $\mathbf{Y}$. The dual distribution of $p(\theta)$, corresponding to the affine subspace $W\hat{=}\mathrm{Ker}\ \{\mathbf{Y}^{\prime\mathrm{T}}\}$, has the analytic form

[TABLE]

where $(\mathbf{s},\Phi)$ are the parameters of the affine subspace $W$.

Proof:

We create an orthogonal basis $\mathbf{u}_{1},\mathbf{u_{2}},\ldots,\mathbf{u}_{M-1}$ for the subspace $V$ and its orthogonal complement $V^{\perp}$ using the Gram-Schmidt orthonormalisation procedure for the projections $\mathbf{p}_{i}=\mathbf{P}\mathbf{e}_{i}$, $i=1,2,\ldots,K$ onto the subspace $V$ (see Appendix B) and similarly for the projections $\mathbf{p}_{i}=(\mathbf{I}-\mathbf{P})\mathbf{e}_{i}$, $i=K+1,K+2,\ldots,M-1$ on its orthogonal complement $V^{\perp}$. We marginalise the $M-1$ dimensional Gaussian over all the parallel affine subspaces of $V$, i.e., those which are of the form $W=\mathbf{w}+V$. We obtain the conditional probability density of $\mathbf{w}\hat{=}\mathbf{s}\in\mathbb{R}^{M-K-1}$

[TABLE]

which is a mean zero, $M-K-1$ dimensional Gaussian with the identity matrix as the covariance matrix.

As to the distribution of directions $\Phi$, isotropic Gaussians imply uniform directions on half hyperspheres. We decompose the vector as $\Phi=(\Phi_{1},\Phi_{2},\ldots,\Phi_{\tilde{K}})$ and further $\Phi_{k}=(\phi_{1}^{k},\phi_{2}^{k},\ldots,\phi_{M-2k}^{k})$. Thus we have

[TABLE]

where

[TABLE]

$k=1,2,\ldots,\tilde{K}$.

So finally,

[TABLE]

and the claim follows. ∎

We have thus derived an analytic form for the probability density of the affine subspace $W=W(\mathbf{s},\Phi)$ or, equivalently, the probability density of the left nullspace of the reduced coefficient matrix $\mathbf{Y}^{\prime}$, assuming the model (16) and Gaussian distributed parameter vector $\theta$.

{unmarginalised}

Let $\mathbf{Z}$ be a $M\times K$ matrix containing linearly independent column vectors $\mathbf{z}_{k}$, $k=1,\ldots,K$, where $K=\mathrm{rank}(\mathbf{Y}^{\prime})$ and $\mathrm{Ker}\ \{\mathbf{Z}^{\mathrm{T}}\}=\mathrm{Ker}\ \{\mathbf{Y}^{\prime\mathrm{T}}\}$. In other words, $W_{0}$ is the intersection of the hyperplanes $\mathbf{z}_{k}^{\mathrm{T}}\theta^{\prime}=\mathbf{0}$, $k=1,2,\ldots,K$ in the dual space, where we parameterise the hyperplanes by their signed distance $\mathbf{s}=(s_{1},s_{2},\ldots,s_{K})$ from the origin and their unit normal vectors $\mathbf{v}_{k}=\mathbf{v}(\varphi^{k})$ using the modified spherical coordinates $\varphi^{k}$, $k=1,2,\ldots,K$, $\Phi=(\varphi^{1},\varphi^{2},\ldots,\varphi^{K})$, as described in the previous section.

The subspace $Z$ spanned by $\{\mathbf{z}_{1},\ldots,\mathbf{z}_{K}\}$ is orthogonal to any affine subspace $W$ parallel to $W_{0}$ in $\mathbb{R}^{M-1}$. These affine subspaces can be hence parameterised by the $K$ dimensional vector $\mathbf{\mathbf{w}}$ representing the point of intersection between $W$ and $Z$. We define $\mathbf{\mathbf{w}}$ in the orthonormal basis $B=\{\mathbf{u}_{1},\ldots,\mathbf{u}_{K}\}$, obtained from the vectors $\mathbf{v}_{1},\ldots,\mathbf{v}_{K}$ by the Gram–Schmidt orthonormalisation procedure. The conditional probability density of $\mathbf{w}$ is obtained as

[TABLE]

which is a mean zero, $K$ dimensional Gaussian with the identity matrix as the covariance matrix. For the Gaussian, we make the coordinate transform $\mathbf{w}=\mathbf{w}(\mathbf{s})$, where

[TABLE]

whereas $\tilde{\mathbf{V}}=(\tilde{\mathbf{v}}_{1}\ \ldots\ \tilde{\mathbf{v}}_{K})$ represent the direction vectors $\mathbf{v}_{1},\ldots,\mathbf{v}_{K}$ in the basis $B$. Hence,

[TABLE]

On the other hand, since the Gaussian distribution is isotropic, the distribution of directions is uniform with the density

[TABLE]

Hence,

[TABLE]

Again, we use the modified spherical coordinates to parameterise the the reduced joint feature vectors, i.e., $\mathbf{Z}=\mathbf{Z}(\mathbf{r},\Phi)$, where $\mathbf{r}=(r_{1},\ldots,r_{K})$ and $\Phi=(\varphi^{1},\ldots,\varphi^{K})$. In addition, we use the fact that $\mathbf{s}=\mathbf{s}(\mathbf{r})$ and $\mathbf{v}_{k}=\mathbf{v}_{k}(\varphi^{k})$, $k=1,\ldots,K$. Now, we get

[TABLE]

where $\mathbf{r}_{-1}=(1/r_{1},1/r_{2},\ldots,1/r_{K})$.

Now, we are interested in the total probability of $\mathrm{Ker}\ \{\mathbf{Y}^{\prime\mathrm{T}}\}$ but there are infinitely many matrices that have the same kernel. Therefore, we need to marginalise over all those $\mathbf{Z}$ that has the kernel equal to $\mathrm{Ker}\ \{\mathbf{Y}^{\prime\mathrm{T}}\}$. Hence,

[TABLE]

where $\mathbf{Z}_{0}$ is any fixed $\mathbf{Z}$.

This is an analytic form for the probability density for the kernel of the reduced joint feature matrix $\mathbf{Y}^{\prime}$ assuming the model (16) and the Gaussian distribution for the parameter vector $\theta$.

5 Mappings to Feature Distributions

In the previous section, we derived the general form for the dual distributions in the function of the parameters of the corresponding affine subspace. Now we discuss how we can extract interesting feature distributions from the dual distributions. The interesting feature distribution often has less parameters than the affine subspace or one may be interested only in certain conditional feature distributions, conditioned on some fixed a subset of the features. In these cases the mapping from the feature distribution to the affine subspace $W$ is not necessarily one-to-one and onto. However, we may form the conditional distribution

[TABLE]

conditioned on the restriction $R\subset\Omega$, where $R$ is parameterised by the interesting part $\tilde{\mathbf{x}}\in\mathbb{R}^{\tilde{N}}$ of the feature distribution.

Theorem 5.1

Let us consider a smooth submanifold $R$ of the parameter space $\Omega$ so that there is a differentiable bijective map between the parameters $\tilde{\mathbf{x}}\in X\subset\mathbf{R}^{\tilde{N}}$ and $(\mathbf{s},\Phi)\in R$ almost everywhere. Then

[TABLE]

where $\mathbf{J}=\frac{\partial(\mathbf{s},\Phi)}{\partial\tilde{\mathbf{x}}}$.

Proof:

Let $\varphi:$ $X\rightarrow R$ so that $(\mathbf{s},\Phi)=\varphi(\tilde{\mathbf{x}})$ is continuous and invertible almost everywhere on $X\subset\mathbb{R}^{\tilde{N}}$. Given an orthonormal basis $B=\{\mathbf{u}_{1},\mathbf{u}_{2},\ldots,\mathbf{u}_{\tilde{N}}\}$ on the tangent space $T_{\varphi_{0}}(R)$

[TABLE]

where $\varphi_{0}=\varphi(\tilde{\mathbf{x}}_{0})$ and $\varphi\in T_{\varphi_{0}}(R)$. The Jacobian of the mapping $\tilde{\mathbf{x}}\mapsto\xi$ is the $\tilde{N}\times\tilde{N}$ matrix

[TABLE]

where we have used the property $\xi_{i}=\mathbf{u}_{i}^{\mathrm{T}}(\varphi-\varphi_{0})$, $i=1,\ldots,\tilde{N}$. Thus

[TABLE]

which holds almost everywhere and is independent of the choice of the orthonormal bases. Using the substitution rule for integrals, the conditional distribution of $\tilde{\mathbf{x}}$ takes the form

[TABLE]

∎

As $p(\tilde{{\mathbf{x}}}|X)$ is easily computed up to a global constant, we may draw samples from it by generating a MCMC chain with the Metropolis–Hastings sampling rule. Moreover, at least in cases where the interesting distribution covers a linear manifold, direct sampling methods can be applied similar to those derived in [12].

6 Experiments

In this section we show two application examples of the dual distributions. We create confidence intervals for conics (Section 6.1) and show how probabilistic point transfer can be constructed by using the covariance information of the estimated trifocal tensor (Section 6.2).

6.1 Dual Distributions for Conics

The points on a conic satisfy the homogeneous quadratic equation

[TABLE]

which is a bilinear equation in $\mathbf{x}$ and $\mathbf{A}$ is a symmetric $3\times 3$ matrix. This equation can be written in the form, using the Veronese embedding,

[TABLE]

where $\theta\hat{=}(a_{11},a_{22},a_{33},2a_{12},2a_{23},2a_{13})$ and $\mathbf{y}\hat{=}(x_{1}^{2},x_{2}^{2},x_{3}^{2},x_{1}x_{2},x_{2}x_{3},x_{1}x_{3})$, i.e., a conic forms a 5-dimensional linear subspace in the six dimensional joint feature space. We assume that the parameter vector estimate $\theta=\theta_{0}$ and its covariance matrix $\mathbf{C}_{\theta}$ are available.

If we dualise the relationship (40) above, we see that a fixed point on the image plane determines a 5-dimensional linear subspace in the dual space, i.e., the space of all conics that intersect the point on the image plane, see Fig. 3. Moreover, to characterise the probability of the point on the image we may construct the total probability (density) of all those conics containing the point by using the uncertainty distribution of the estimated conic. In other words, the dual distribution characterises the confidence of the estimated conic by illustrating what has been learned from the locations of the points on the true conic.

By the way of an example, we estimated the maximum likelihood estimate and its covariance matrix for the conic containing 25 points, shown in Fig. 4a, assuming i.i.d. Gaussian noise in the 2D measurements. To evaluate the dual density at the selected location $(x_{1},x_{2})$ on the 2D plane, we need to evaluate (15) as well as compute the right magnification factor. As we parameterised the reduced joint feature vector $\mathbf{y}^{\prime}$ by the modified spherical coordinates, we may construct the mapping $(x_{1},x_{2})\mapsto(\rho,\phi)$ and its Jacobian to evaluate the dual pdf using (38). The dual pdf for the estimated conic is illustrated in Fig. 4b. The contours visualise the fact that the we have a strong belief about the true conic points near the training data but extrapolation beyond the training data contains a substantial risk.

{conic}

Using the asymptotic normality property of the MLE, we assume that $\theta\sim N(\theta_{0},\mathbf{C}_{\theta})$ with the density function $p(\theta)$. Since $p(\theta)$ is Gaussian, we may analytically marginalise it over the hyperplanes $S_{\mathbf{n}}$ parallel to $\mathbf{y}$ in the dual space, i.e.,

[TABLE]

where $\mathbf{n}$ denotes the unit normal of the parallel hyperplanes and $\lambda$ parameterises the (negative) signed distance of a parallel hyperplane from the origin in the dual space, thus we define $\mathbf{y}=(\mathbf{n}^{\mathrm{T}}\ \lambda)^{\mathrm{T}}$. Since we marginalise the 6-dimensional Gaussian distribution over the 5-dimensional subspace, $p(\lambda|\mathbf{n})$ is one-dimensional normal distribution.

We may also derive the distribution of directions that define the orientation of $\mathbf{n}$. This is equivalent in marginalising the 6-dimensional Gaussian distribution over the signed radius in the modified111Instead of defining the radius parameter on the half infinite interval $[0,\infty)$, we define it on the full interval $(-\infty,\infty)$, and redefine the angle domains accordingly. 6-dimensional spherical coordinates as

[TABLE]

where $\Phi$ denotes variables in the modified 6-dimensional spherical coordinate system. This integral can be analytically computed, as will be seen in next section.

The two distributions above can be combined to form the distribution for $\mathbf{y}$

[TABLE]

The remaining problem is to transform the density of $\mathbf{y}$, parameterised by $\mathbf{n}$ and $\lambda$, onto the inhomogeneous image coordinates $(u,v)$. The difficulty is that the joint feature mapping $\mathbf{y}=\mathbf{y}(\mathbf{x})$ is not one-to-one and onto. Therefore every $\mathbf{y}$ does not correspond to a true image point. However, since the uncertainty of the the conic reflect the behaviour of the training data, the variation of the conic parameters should be small in the directions of forbidden $\mathbf{y}$, so those can be neglected. We thus make an intermediate transformation $\mathbf{z}(\mathbf{y})=(y_{1},y_{2},y_{3},y_{4},y_{4}/y_{5},y_{4}/y_{6})$, hence, $(u,v)\equiv(z_{5},z_{6})$. Now, we may write

[TABLE]

where we have marginalised the “nuisance” parameters $z_{1},\ldots,z_{4}$ out from the distribution.

Sampling from the dual distribution (44) is straightforward. First we draw a sample from the six-dimensional Gaussian $p(\theta)$ and compute the direction vector $\mathbf{n}$ of the sample. Given $\mathbf{n}$, we then draw a sample for $\lambda$ from the 1D Gaussian (41). Since, $\mathbf{y}=(\mathbf{n},\lambda)$, we finally obtain the estimate for the point $(u,v)$, as $(y_{4}/y_{5},y_{4}/y_{6})$.

6.2 Probabilistic Point Transfer with an Uncertain Trifocal Tensor

The geometry of three projective views is characterised by several incidence relations or trilinearities, which are collected into Table I [1]. According to our preferences, we could use any of the trilinearities to create dual distributions or conditional dual distributions. As an example, we now illustrate how the trifocal point transfer can be “probabilised” by constructing a dual distribution for the trifocal tensor given estimates for the tensor and its covariance matrix as well as a novel point match in two views. In this case, the dual distribution is the conditional position distribution of the transferred point in the third view. We use the nine point–point–point constraint equations (the first relation in Table I) of which only four equations are independent. Therefore the dimension of the affine subspace of this example is $K=(M-1)-4=14$.

Given the point match $\mathbf{m}^{\prime}\leftrightarrow\mathbf{m}^{\prime\prime}$ in the views two and three, the construction of the conditional probability density $p(\mathbf{m}|\mathbf{m}^{\prime},\mathbf{m}^{\prime\prime},\mathcal{T},\mathbf{C}_{\mathcal{T}})$ is as follows. The trilinear relations for three points define the elements of $\mathbf{Y}$ in the model (16). After constructing $\mathbf{Y}$, the whitened matrix $\mathbf{Y}^{\prime}$ is obtained from (17). However, we regularised the whitening transform by replacing $\tilde{\Lambda}$ by $\tilde{\Lambda}+\lambda\mathbf{I}$, where $\lambda=10^{-8}$ since the number of data points was limited in our experiment and hence, due to overfitting of the covariance matrix, the smallest ones of the eighteen non-zero eigenvalues were (assuming a non-degenerate configuration) relatively small.222The diagonal Tikhonov regulariser was adequate here but a more elegant way would have been investigating the numerical rank of the covariance matrix and developing an automatic model selection scheme to determine the dimension of the affine subspace. Degenerate configurations could be also handled in this way. The probability density for the affine subspace, corresponding to $\mathbf{Y}^{\prime}$, is given by (25) and the desired conditional probability density values are finally obtained from (38), up to scale. In short, to compute the conditional probability density value for a location of interest in the first view, the corresponding affine subspace parameters are computed from $\mathbf{Y}^{\prime}$ and the probability density value of the corresponding affine subspace is weighted by the term containing the Jacobian of the transformation.

To illustrate the conditional point transfer density, we estimated the maximum likelihood trifocal tensor and its covariance matrix from 274 outlier free point correspondences for the image triplet shown in Fig. 5 assuming i.i.d. Gaussian noise in the measurements [1]. Then, by using a novel point correspondence in the views two (Fig. 6b) and three (Fig. 6c), we computed the conditional probability density in the first view, as reported above. To visualise the shape of the pdf, we selected, three pdf values at the levels $10^{-1}$, $10^{-2}$, and $10^{-3}$ times the maximum value and show the corresponding contours in Fig. 6a and 6d. It can be seen that the transferred density is non-Gaussian while it indicates the feasible locations for the correspondence. The pdf has its maximum close to the point where the ML epipolar lines meet whereas the local shape of the peak is oriented towards the mean axis of the two epipolar lines. The pdf shape also seems to illustrate the well known fact that the trifocal constraint is more versatile than the mere epipolar geometries between the three views.

7 Discussion

The dual distributions are a tool for Bayesian inference with uncertain multilinear models, used e.g. in geometric image analysis. According to the Bayesian paradigm, we handle geometric relationships as probability distributions in contrast to the traditional approaches, which are based on deterministic estimates. Moreover, the theory provides tools for propagating the geometric uncertainty information to the dual variables. This approach helps especially in the cases when there is no explicit form for the mapping but a multilinear incidence relation which indirectly connects the parameter space and the observations.

As in the point–line case there are certain points and lines with a special role [12], likewise there are special subspaces in the general case to be investigated in future. They arise from the affine subspaces spanned by the eigenvectors of the covariance matrix—the deterministic trifocal point transfer, for instance, should be constructed by identifying the most likely affine subspace in the Gaussian model using the point–point–point incidence relation. In addition, there are many other ways to utilise the special subspaces arising from the multiple constraint equations. However, one should be careful in the interpretations as covariance estimates are not generally invariant to the selection of the coordinate frame.

This work is fundamentally a study of the generalised Radon transform for probability measures. Since the transform is based on integration over all the affine subspaces, a suitable parameterisation for affine subspaces is required. A natural algebraic choice would be the Grassmannian parameterisation, and it is likely that it would provide certain algebraic benefits that are to be investigated in future. For instance, one should study the transforms with respect to a natural Riemannian metric, invariant to motions of the coordinate system, to make the approach invariant to Gauge transforms. Proceeding in this direction would clarify the role of the dual distributions in between Santalo’s and Gelfand’s schools of integral geometry.

As far as the numerical development is concerned, more efficient numerical tools would be advantageous for geometric inference. The study of special subspaces and efficient sampling from the dual distributions in the general case are additionally to be investigated in future. Likewise, the study of marginal distributions, and proper handling of nonlinear geometric constraint manifolds are part of the future work.

8 Conclusions

In this paper we have shown how the parameter distributions of multilinear models can be dualised. The proposed approach provides means for a pure statistical treatment of multilinear relations where the uncertainty of the geometric entity is taken into account. The dual distributions are closely connected to integral geometry and have analytic forms with relatively small assumptions. We demonstrated the applicability of the theory by characterising the confidence of estimated conics and constructed the probabilistic trifocal point transfer by using an uncertain trifocal tensor with its estimated covariance matrix. In future, an interesting research direction is to investigate the dual distributions from the view point of integral geometry to deeply understand the theoretical connections as well as to develop additional numerical methods for applications.

A. Modified Spherical Coordinates in $\mathbb{R}^{N}$

In this paper, we use modified spherical coordinates in $N$-dimensional space. That is, we assign a sign for the radius in the conventional $N$-dimensional spherical coordinates [19] and parameterise the directions using points only on the other half of the unit hypersphere. This parameterisation is natural for one-dimensional subspaces in $\mathbb{R}^{N}$ as the directions then uniquely parameterise the linear subspaces almost everywhere and the (signed) radial parameter parameterises the points in the subspace.

Assume that333If $N=1$ or $N=2$ the construction is similar but these cases are omitted here for simplicity. $N\geq 3$ and let $\rho\in\mathbb{R}$ be the radial coordinate and $\phi_{1},\phi_{2},\ldots,\phi_{N-1}$ be the angular coordinates so that $\phi_{1}$ takes values between 0 and $\pi/2$, $\phi_{2},\ldots,\phi_{N-2}$ are between [math] and $\pi$ and $\phi_{N-1}$ is between [math] and $2\pi$. The Cartesian coordinates $x_{i}$, $i=1,\ldots,N$ are then defined as

[TABLE]

The inverse transform is

[TABLE]

and the volume element is obtained as

[TABLE]

B. Proof of Lemma 1

Proof:

It is well known that an orthogonal projection matrix $\mathbf{P}$ may be decomposed as $\mathbf{P}=\mathbf{U}\mathbf{U}^{\mathrm{T}}$ where the columns of $\mathbf{U}$ form an orthogonal basis of the range of $\mathbf{P}$. We hence form an orthonormal basis for the range from the basis $\mathbf{p}_{i}=\mathbf{P}\mathbf{e}_{i}$, $i=1,\ldots,K$ by using the Gram–Schmidt orthonormalisation procedure. Now, let

[TABLE]

and let $\mathbf{L}=(\mathbf{u}_{1}\ \mathbf{u}_{2}\ \ldots\ \mathbf{u}_{K})$, i.e., the columns of $\mathbf{L}$ form an orthonormal basis for the range of $\mathbf{P}$, and $\mathbf{P}=\mathbf{L}\mathbf{L}^{\mathrm{T}}$.

We need to show that $\mathbf{L}$ is additionally a lower triangular matrix which has a strictly positive diagonal. We note that

[TABLE]

where the second term in the nominator is zero because the vectors $\mathbf{u}_{k}$, $k=1,2,\ldots,K$ are linearly independent by assumption. It follows that $u_{1k}=0$, $k=2,\ldots,K$. Thus, $u_{11}/|u_{11}|=1\implies u_{11}>0$. Similarly,

[TABLE]

that implies $u_{2k}=0$, $k=3,\ldots,K$ and $u_{22}>0$. Likewise, $u_{lk}=0$, and $u_{ll}>0$ when $k=l+1,\ldots,K$, hence, $\mathbf{L}$ is a lower triangular matrix with strictly positive diagonal.

We still need to show that the representation is unique. Let us assume the contrary, i.e., there are two different lower triangular matrices $\mathbf{L}_{1}=(\mathbf{u}_{1}\ \mathbf{u}_{2}\ \cdots\ \mathbf{u}_{K})$ and $\mathbf{L}_{2}=(\mathbf{v}_{1}\ \mathbf{v}_{2}\ \cdots\ \mathbf{v}_{K})$ which have the the indicated properties. Then $\mathbf{P}=\mathbf{L}_{1}\mathbf{L}_{1}^{\mathrm{T}}=\mathbf{L}_{2}\mathbf{L}_{2}^{\mathrm{T}}$ and

[TABLE]

Since $\mathbf{u}_{1}$ and $\mathbf{v}_{1}$ are unit vectors, it follows that $|u_{11}|=|v_{11}|$. The diagonal elements are additionally strictly positive so $u_{11}=v_{11}$ that implies $\mathbf{u}_{1}=\mathbf{v}_{1}$. Similarly, we see that $\mathbf{u}_{l}=\mathbf{v}_{l}$, when $l=2,\ldots,K$, i.e., $\mathbf{L}_{1}\equiv\mathbf{L}_{2}$ that contradicts the assumption, hence, $\mathbf{L}$ must be unique. ∎