Problem regarding compiling a program- Stiefel Manifold

[rajatsanya...@gmail.com]

First of all, I would like to convey my gratitude to the MANOPT's authors. MANOPT is a great initiative, and makes ours work a lot easier.

I am trying to solve a problem on Stiefel manifold using MANOPT but facing some errors. I could not figure out how solve the problems.

Any generous help or suggestion is highly appreciated. The problem that I want to implement, is as follows-

\begin{equation*}
\min_{O_1, \ ...\ , \ O_M \in \mathbb{O}(d)} \Tr (C O^T O)
\end{equation*}
\begin{equation*}
\text{where, } O=[O_1, \ ...\ , \ O_M ], \ C \in \mathbb{S}^n_+.
\end{equation*}

A part of my program-

%% Manopt Code ---->
% Create the problem structure.
manifold = stiefelfactory(d,d,M);
problem.M = manifold;
% Define the problem cost function and its Riemannian gradient.
problem.cost = @(O) cost(O);
function f = cost(O)
f = 0;
for i = 1:M
for j = 1:M
f = f+trace(O(:,:,i)*C((i-1)*d+(1:d),(j-1)*d+(1:d))*...
O(:,:,j)');
end
end
end
problem.grad = @(O) manifold.egrad2rgrad(O,egrad(O));
function g = egrad(O)
g = zeros(d,M*d);
for i = 1:M
h = zeros(d,d);
for j = 1:M
h = h+2*O(:,:,j)*C((j-1)*d+(1:d),(i-1)*d+(1:d));
end
g(:,(i-1)*d+(1:d)) = h;
end
end
% Solve.
[O_s,f_cost,info,options] = trustregions(problem);

[BM]

Dear Rajat,

Thank you for your comments.

I am a bit confused by the cost function that you provided. I do not see O_1, O_2 in the cost function in the following problem formulation.  I shall be glad if you could be precise here as this will allow members to give concrete feedback.

\begin{equation*}
\min_{O_1, \ ...\ , \ O_M \in \mathbb{O}(d)} \Tr (C O^T O)
\end{equation*}
\begin{equation*}
\text{where, } O=[O_1, \ ...\ , \ O_M ], \ C \in \mathbb{S}^n_+.
\end{equation*}

Also, what is \mathbb{S}^n_+? If these are rotations, it might be better to use the rotations factory.  

Regards,

BM

[Nicolas Boumal]

Dear Rajat,

Thank you for your message. As I understand it, the problem you want to solve is what people sometimes call orthogonal synchronization (I imagine you know this paper, but here is a nice convex relaxation for it: http://arxiv.org/abs/1308.5207.)

I happen to have considered that problem quite a bit myself. I invite you to look at this paper:

A Riemannian low-rank method for optimization over semidefinite matrices with block-diagonal constraints

http://arxiv.org/abs/1506.00575

Code (using Manopt) is available here:

https://perso.uclouvain.be/nicolas.boumal/Staircase/RiemannianStaircase_v1.2.zip

For your usage, I would recommend you take a look at the file stiefelstackedfactory in that zip file: it will allow you to optimize over St(d, d)^M in the form that you want to consider (i.e., with the O_i's stacked on top of each other). That factory will be part of the next release of Manopt (very soon now!).

Note that this code proposes to optimize not over O(d)^M, but rather over Stiefel(d, d')^M, with d' > d. This is because O(d) has two connected components, so that O(d)^M has 2^M components, and it's almost impossible to know which one to explore. (As the paper explains, very often, we can get back the O(d) solution in the end.) A particular case where we do know, is when we actually want to optimize over SO(d)^M, that is, orthogonal matrices with determinant 1. In that case, you may use the rotationsfactory (geometry for SO(d)) in Manopt. I have a short paper about that here:

Robust estimation of rotations from relative measurements by maximum likelihood

https://perso.uclouvain.be/nicolas.boumal/papers/Boumal_Singer_Absil_Robust_estimation_of_rotations_from_relative_measurements_by_maximum_likelihood_2013.htm

with code (using Manopt) here (it uses an older version of Manopt, but should be plug and play with the newer versions)

https://perso.uclouvain.be/nicolas.boumal/SynchronizeMLE/SynchronizeMLE_1.zip

Bamdev: \mathbb{S}^n_+ is the set of positive semidefinite matrices. It really doesn't matter though, since the cost is invariant under diagonal shifting of C, so even if C is not psd, we can always shift it to make it psd, by adding lambda_min * identity. The distinction mostly matters for researchers establishing approximation bounds, where the actual value of the cost is important, not just the arg-min or arg-max.

If you would like more help with your code, could you please give a complete test file so that we can see what the errors seems to be? Without a description of the errors / failures, it's hard to tell.

Best,

Nicolas

[rajatsanya...@gmail.com]

Thank you for your quick response.
I want to mean that C is a fixed matrix which belongs to the set of PSD matrices.
I want to create a manifold by using-

manifold = stiefelfactory(d,d,M),

but in the documentation I noticed that-

% Points are represented as matrices X of size n x p x k (or n x p if k=1,
% which is the default) such that each n x p matrix is orthonormal,
% i.e., X'*X = eye(p) if k = 1, or X(:, :, i)' * X(:, :, i) = eye(p) for
% i = 1 : k if k > 1. Tangent vectors are represented as matrices the same
% size as points.

So, I thought that O_i could be presented as O(:,:,i). Please correct me if I was wrong. Could you please give me some suggestions about how to write the desired cost function that? I think that I might do the same mistake while computing the gradient.

[rajatsanya...@gmail.com]

Dear Prof. Boumel,

Thank you sir for your insight. I got the point. I will definitely go through the material that you have referred.

Thank you again for your generous help.

Sincerely,
Rajat.