New Manifold: Symmetric Unitary Matrices

[Leonardo Mörlein]

Hi,

first of all, I want to thank you for the great tool that you have developed with manopt. I really like the way it is designed! :)

Now coming to the reason for my post: I currently want to optimize something on a manifold that is not yet part of manopt. It is the manifold of Symmetric Unitary Matrices, which occurs e.g. in circuit/scattering theory and represents circuits that are both lossless and reciprocal.

Doesn't really matter from a mathematical point of view, but for the background:

Condition for losslessness is:

Condition for reciprocity (which almost all circuits should fulfill) is:


Now, I want to add this new manifold and I am still struggeling a bit how to do this. My current approach was to try to (naively) copy most of the stuff from the unitaryfactory.

My current tests look somehow like this:

gsmfact = unitaryfactory(N+P, K);
gsmfact.retr = @(Theta, U, t) multisym(gsmfact.retr(Theta, U, t));
gsmfact.proj = @(X, H) multisym(gsmfact.proj(X, H));
gsmfact.egrad2rgrad = gsmfact.proj;
gsmfact.rand = @() multisym(gsmfact.rand());

However, it doesn't seem to work as I expect. Probably, I missed to adjust some of the fields or adjusted them incorrectly. The optimization stops after two steps due to small step size.

About my background: I am electrical engineer, not very deep into the theory behind manifolds, but have a rough idea what happens and what the terms mean.

Now, coming to my questions:

1) Is there a guide that I can follow how to add new manifolds?

2) Can you maybe point me towards which "methods" of the unitaryfactory() manifold I have to adjust and which I ones can leave as they are?

2a) I think, I can keep dim(), inner(), norm(), typicaldist(),  egrad2rgrad(), hash(), lincomb(), zerovec(), transp(), isotransp(), pairmean(), dist(), vecmatareisometries() and lie_identity(). Do you think, this could be correct?

2b) I think, I have to change proj(), rand(), randvec() and retr(). Do you think, this could be correct?

2c) I am currently unsure about: tangent(), tangent2ambient(), ehess2rhess(), exp(), log(). What do you think? Do I have to change them?

3) I saw the remark that tangent vectors are represented as skew-hermitian matrices in the unitaryfactory and tangent2ambient() is necessary to switch to the embedding space.

3a) I guess, this means that simply wrapping the .proj() function in multisym() is not a good idea, since "symmetrizing" works differently in the Lie algebra than in the embedding space, doesn't it?

3b) Do I interpret it correctly that M.tangent() is the inverse operation of M.tangent2ambient()? If so, I could use this to check 3a), couldn't I?

3c) Can I keep the representation in that Lie algebra? Or is another representation/another Lie algebra necessary for symmetric unitary matrices than for unitary matrices in general?

3b) If it is not too complicated, can you maybe explain in a few simple words why/if the representation in the Lie algebra is necessary? Or point me to a reference where I can understand better what is going on here?

4) Would the manifold of unitary symmetric matrices be of general interest to be added to the manopt code?

Not sure if it is helpful/necessary, but if you wish, I could also think about a minimal cost function that could be used with this manifold. (My current cost function is rather complex and uses a product of different manifolds.)

During writing this post, I just found checkmanifold(). In the meantime while waiting for an answer from you, I will try to come forward using checkmanifold(). I will let you know about my progress.

Also: If you think, the idea that I could implement this myself without a lot of deep background knowledge about manifolds and Lie Algebras is a little too ambitious, feel free to let me know.

Cheers,
Leo

[Leonardo Mörlein]

Well, I just noticed that calling multisym() on an unitary matrix makes the result non-unitary. I will have to look for another way.

Cheers,
Leo

[Leonardo Mörlein]

I got something to "work":

gsmfact = unitaryfactory(N+P, K);

function Y = multipoldecomp(X)
    for kk = 1:size(X, 3)
        Y(:,:,kk) = poldecomp(X(:,:,kk));
    end
end

symunit = @(x) multipoldecomp(multisym(x)); % takes unitary x as input and symmetrizes it.

gsmfact.retr = @(Theta, U, t) symunit(gsmfact.retr(Theta, U, t));
gsmfact.rand = @() symunit(gsmfact.rand());

function U = randomvec(gsmfact, x)
    symunit = @(x) multipoldecomp(multisym(x)); % takes unitary x as input and symmetrizes it.

    randvecunitary = gsmfact.randvec(x);
    randvecunitaryambient = gsmfact.tangent2ambient(x, randvecunitary);
    randvecunitaryambientsym = symunit(randvecunitaryambient);
    U = gsmfact.proj(x, randvecunitaryambientsym);
    nrmU = sqrt(U(:)'*U(:));
    U = U / nrmU;
end
gsmfact.randvec = @(x) randomvec(gsmfact, x);

I took code for polar decomposition (poldecomp) from here: https://de.mathworks.com/matlabcentral/fileexchange/48735-polar-decomposition

The idea is that I use multisym() to make the unitary matrix symmetrical and poldecomp makes it unitary (again).

This seems to work (somehow). checkmanifold() is happy. And the result vectors are indeed symmetric unitary.

However, I am not sure whether y = symunit(x) returns the y, which is the closest possible symmetric unitary matrix to the unitary matrix x. Probably not? And is this even important? I don't know...

Cheers,
Leo

[Nicolas Boumal]

Hello Leo,

This is an interesting set -- had not seen this one before.

Just to make sure I got this correctly: (in Matlab notation), we're considering the set of complex matrices S of size n-by-n such that S'*S = eye(n) (unitary)  and  S.' = S (so, both real(S) and imag(S) are symmetric).

Do we know for sure this is a smooth manifold, embedded in C^{n \times n}?

If I split S = A + i B where A, B are real matrices of size n-by-n, then S'*S = A'*A + B'*A + i ( A'*B - B'*A ), which is identity (that is, S is unitary) iff A'*A + B'*B = eye(n), and A'*B = B'*A.

Then also requiring S.' = S effectively forces A = A' and B = B' (both real symmetric), so the complete list of requirements becomes:

• A = A'
• B = B'
• A'*A + B'*B = eye(n)
• AB = BA

The last one says that A and B commute. Since they are symmetric matrices, this means that A and B can be diagonalized by a common orthogonal matrix of eigenvectors: let's call it V.

Thus, A = VDV' and B = VEV', where D and E are two real, diagonal matrices, and V is a real, orthogonal matrix. The only other condition is that D² + E² = I, which effectively means that

D = diag(cos(theta_1), ..., cos(theta_n)) and E = diag(sin(theta_1), ..., sin(theta_n)).

Overall, it seems to me that your set can be smoothly parameterized by taking n points on a circle (the angles theta_1, ..., theta_n) and a real, orthogonal matrix V, then forming A and B as above and assembling them to form S.

In Manopt, this might go as follows, using a product manifold:

elems.angles = obliquefactory(2, n);  % could also use complexcirclefactory(n) if you're comfortable with mixing complex numbers in there

elems.basis = rotationsfactory(n);     % could also write stiefelfactory(n, n) if you run into issues related to representation of tangent vectors

manifold = productmanifold(elems);

X = manifold.rand();

S = X.basis*diag(X.angles(1, :) + 1i*X.angles(2, :))*X.basis';  % if using complexcirclefactory, here we'd just take diag()

This indeed generates a random matrix S, and we have a smooth parameterization for that set. You could write a cost function f(S) as a function g(angles, basis) = f(S) (with the change of variables shown above), and get the graident/Hessian using manoptAD (auto-diff) or working it out by hand. This would work.

Or you could implement a factory for that set of S directly, as you started doing above. There, it's indeed important to makes sure the retraction is correct. Even before doing so, one would want to make sure that the set of S is indeed smooth (the above is a smooth parameterization, but that does not imply the set itself is smooth). It very well may be, but it's not immediately obvious to me.

Best,

Nicolas

[Nicolas Boumal]

The version with complexcirclefactory and stiefelfactory (you can mix and match) looks like this:

elems.angles = complexcirclefactory(n);  % vector of n unit-modulus complex numbers
elems.basis = stiefelfactory(n, n);      % manifold of n-by-n real, orthogonal matrices
manifold = productmanifold(elems);

If you test it out:

X = manifold.rand();
S = X.basis * diag(X.angles) * X.basis';

norm(S'*S - eye(n))
ans =
   1.7494e-15
norm(S.' - S)
ans =
   1.8618e-16