How to solve least square problems

[Markus Sprecher]

Is there a nice way to solve least square problems with manifold-valued data? Second, it also possible to estimate the uncertainties of the reconstructed data?

Consider for example the following model problem: We are given n points u_1,...,u_n in R3 and measurements v_1,...,v_n in R3 of the image of these points under an orientation preserving isometry. The goal is find the isometry (rotation and translation). I was able to write a code (see below) which solves this model problem. However, I am struggling when the problem becomes more complicated. I am also not able to provide the Hessian of my costfunction. Is there maybe a possibility to provide manopt the residual and manopt can compute the derivatives of the cost function (=Frobenius norm of residual)?

Furthermore I would also like to provide some uncertainty measurements. If I consider my measurements V as a random Variable I can consider the reconstruction also as a random variable. I want to estimate the variance of the reconstruction given the variance of V.


function [I]=manopt(U,V)

%number of points
n=size(U,2);
problem.M=specialeuclideanfactory(3);

%residual
res=@(X) X.R*U+X.t*ones(1,n)-V;

%costfunction
problem.cost=@(X) .5*sum(sum(res(X).^2));
%Gradient
problem.egrad=@egrad
function G=egrad(X)
   G.R=res(X)*U';
   G.t=sum(res(X),2);
end

%checkgradient(problem)

% Solve.
I = trustregions(problem);

end

[Nicolas Boumal]

Hello Markus,

Thank you for your message.

 I am struggling when the problem becomes more complicated.

What kind of complications do you envision? Going to a different set of transformations, or optimizing a different cost function?

 

I am also not able to provide the Hessian of my cost function.

As is, Manopt will approximate the Hessian with finite differences, automatically, if the Hessian is not provided by the user. Typically, this works quite nicely. In your code, when you call trustregions, this is what happens under the hood.

This being said, if you want to provide the Hessian in closed form as well, the easiest way is to provide the Euclidean Hessian. Then, problem.ehess = @ehess, where ehess is a function of 2 inputs: X, a structure containing R and t, is the root point; and Xdot is a structure with also fields R and t, which are the directions in which X.R moves and X.t move. To obtain expressions for the output, say H, also a structure with fields R and t, you need to compute the directional derivatives of G.R and G.t (the gradient) when moving in the direction Xdot.R and Xdot.t. See http://www.manopt.org/manifold_documentation_rotations.html for how to deal with the rotations, in particular.

(There is some subtlety with how tangent vectors are represented, as skew-symmetric matrices. In particular, you probably want to start the code in ehess with "Xdot = problem.M.tangent2ambient(X, Xdot);")

The code is below.

 

Is there maybe a possibility to provide manopt the residual and manopt can compute the derivatives of the cost function (=Frobenius norm of residual)?

That is not included, unfortunately. This is part of a bigger missing feature: to have automatic differentiation.

I would also like to provide some uncertainty measurements. If I consider my measurements V as a random variable, I can consider the reconstruction also as a random variable. I want to estimate the variance of the reconstruction given the variance of V.

This exceeds Manopt's purposes, as it becomes quite specific to the particular problem at hand. One simple way of estimating the variance though is as follows:

 * for repeated realizations of the random noise variable

 * compute the estimator X_estim

 * compute the error vector: E = problem.M.log(X_true, X_estim)

Average the error vectors E over the repetitions to estimate the bias (a tangent vector at X_true) ; remove the bias, then average the norms of estimate the variance. More about bias and variance on manifolds can be found here:

1) http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=6418045 

2) http://imaiai.oxfordjournals.org/content/early/2013/09/23/imaiai.iat006.abstract  also here:  http://arxiv.org/abs/1211.1621 

(Ref. 2 focuses on rotations.)

 

Best,

Nicolas

function [I]=foo(U,V)

n=size(U,2);

problem.M=specialeuclideanfactory(3);

% Residual

res = @(X) X.R*U + X.t*ones(1,n)-V; 

% Cost function

problem.cost=@(X) .5*norm(res(X), 'fro')^2;

% Gradient

problem.egrad = @egrad;

function G=egrad(X)

   r = res(X);

   G.R = r*U';

   G.t = sum(r, 2);

end

% Hessian

problem.ehess = @ehess;

function H = ehess(X, Xdot)

   Xdot = problem.M.tangent2ambient(X, Xdot);

   dres = Xdot.R*U + Xdot.t*ones(1,n);

   H.R = dres*U';

   H.t = sum(dres, 2);

end

% checkgradient(problem); pause;

% checkhessian(problem); pause;