🚀 go-pugleaf

I would appreciate greatly any help in understanding the following
matrix factorization problems.

If matrices A and B are real and skew-symmetric (e.g.
A=-A(transpose)), then the matrices exp(A) and exp(B) are orthogonal
with determinants of +1.  The product U=exp(A)*exp(B) is also
orthogonal with det(U)=+1, and there exists a real skew-symmetric
matrix C which satisfies exp(C)=U.  The matrix C is not unique, but
this is OK.

The cases I am interested in are such that A and B do not commute and
have a kind of "complementary" structure: If A(i,j) is nonzero, then
B(i,j)=0 for row and column indices "i" and "j"; similarly, if B(i,j)
is nonzero, then it follows that A(i,j)=0.  The matrix dimensions of
interest are in the range of 50 to 500 (i.e. N**3 type algorithms are
alright provided the leading coefficient is relatively small, <10)

One special case of this "complementary" structure is when A assumes a
block-off-diagonal form, and B assumes a block-diagonal form.  This is
special because matrices of the form exp(B) form a group and are
closed under multiplication.  By this I mean the product
exp(B)*exp(B') may always be written as exp(B") where all three
matrices B, B', and B" have the same block-diagonal form.

I am currently using the routines trizd(), imzd(), and tbakzd() from
TOMS algorithm 530 to perform the transformation A --> exp(A).  This
involves real arithmetic only.  I am satisfied with the performance
and accuracy of this approach, but would welcome any suggestions for
improvements.

********************

Problem 1: What is the best way given an orthogonal matrix U to
determine the matrix C that satisfies U=exp(C)?  Is it possible to do
this with real arithmetic?

[comments: I know how to do this with complex arithmetic by
diagonalizing U and computing X*ln(D)*X(transpose).  The eigenvectors
X and the eigenvalues D are all complex.  Since the reverse
transformation can be done using real arithmetic only, and since the
final result is real, to within machine precision, I suspect there is
a cheaper way of doing this which involves real arithmetic only.  I
realize that C is not unique; I need only one of the infinite number
of solutions.]

********************

Problem 2: Given the matrix U and a "structure" for the matrices A and
B (i.e. the locations of the potentially nonzero elements), what is
the best way to compute A and B which satisfy the relation
U=exp(A)*exp(B)
    a) in the general case where A and B have "complementary" structures?
    b) in the case that A is block-off-diagonal and B is block-diagonal?

[comments: Again, I realize that A and B are not unique; I need only
to find one of the infinite number of solutions.  I have an iterative
scheme to solve case b) above, but it can be slowly convergent and I
expect that for particular choices for U it will not be practical.  This
scheme uses the closure property of the B factors, and therefore does
not generalize to the general case a).  However, I would be interested
in other schemes that also do not generalize.]

********************

Out-and-out solutions to these problems, pointers to software,
literature references, or other constructive comments would be
appreciated.

Ron Shepard