Computing an orthogonal matrix subject to linear constraints

I am looking for a method to solve the matrix equation $$DXa = Xb$$ where $D\in \mathbb{R}^{n\times n}$ is diagonal, $a, b\in \mathbb{R}^{n}$ and $X$ is the unknown orthogonal $n\times n$ matrix such that $X^TX = I$. The matrices $a$ and $b$ have unit norm, and $D$ has diagonal entries of unit magnitude, so that the matrix equation is consistent.

• This looks heavily underdetermined, as you want to solve $n$ equations for $n(n-1)/2$ variables. Is what you wrote really what you want? Nov 19 '12 at 11:57
• Yes, it is underdetermined. For now, I am interested in finding any solution. Nov 19 '12 at 13:09
• Why do your conditions on $a$, $b$, and $D$ imply the consistency of the equations? For $n=1$ consistency already requires the additional condition $Da=b$. Nov 19 '12 at 15:14
• All I meant was I believe that there should always be a solution. For me D is in practice something like +1, -1, +1, -1, ... on the diagonal. Nov 19 '12 at 20:32

Generalizing Arnold's example, suppose $D = I$. Then the problem is to find an orthogonal $X$ satisfying $X(a-b) = 0$, which is possible only when $a = b$.

Interestingly, in 2D with $D = \text{diag}(1, -1)$, one can always find a plane rotation $X$ that does the job. Evidently the solvability for arbitrary unit-norm $a$, $b$ depends only on $D$. What is the general condition that makes the problem solvable? Is it something like $D$ has negative determinant, i.e., it has an odd number of -1's on the diagonal?

• In your example, if we generalize from dimension two to dimension $2n$, with $D = \mathrm{diag}(I,-I)$, is it still always possible? Because that is essentially the problem I have (I think for me the dimension is actually $2n+1$ with the $I$ block being of size 1 larger than the $-I$ block). Nov 22 '12 at 7:36

We have a system of two equations in $$\mathrm X \in \mathbb R^{n \times n}$$

\begin{aligned}\rm D X a &= \rm X b\\ \rm X^\top X &= \rm I_n\end{aligned}

The convex hull of the orthogonal group $$\mathrm O (n)$$ is defined by $$\mathrm X^\top \mathrm X \preceq \mathrm I_n$$, or, equivalently, by the inequality $$\| \mathrm X \|_2 \leq 1$$. Hence, we obtain the following (convex) feasibility problem

$$\begin{array}{ll} \text{minimize} & 0\\ \text{subject to} & \rm D X a = \rm X b\\ & \| \mathrm X \|_2 \leq 1\end{array}$$

Let us try to minimize the trace of $$\rm X$$ instead.

$$\begin{array}{ll} \text{minimize} & \mbox{tr} (\mathrm X)\\ \text{subject to} & \rm D X a = \rm X b\\ & \| \mathrm X \|_2 \leq 1\end{array}$$

Numerical experiment

Using NumPy to randomly generate $$\rm a$$, $$\rm b$$ and $$\rm D$$ and CVXPY to solve the convex program,

from cvxpy import *
import numpy as np

np.set_printoptions(linewidth=250)

n = 16

a = np.random.rand(n)
b = np.random.rand(n)
D = np.diag(2 * np.random.randint(2,size=n) - 1)
I = np.identity(n)

X = Variable(n,n)

# define optimization problem
prob = Problem( Minimize(trace(X)), [ D*X*a == X*b, norm(X,2) <= 1 ] )

# solve optimization problem
print prob.solve()
print prob.status

# print results
print "X = \n", X.value
print "Spectral norm of X = ", np.linalg.norm(X.value,2)
print "Round(X.T * X) = \n", np.round((X.value).T * (X.value), 1)

which outputs the following

-14.9470050407
optimal
X =
[[ -8.90251732e-01   2.33984981e-01   1.70258281e-01  -5.53391898e-02   2.34985256e-02   6.37982593e-02   2.18083930e-01   2.82596557e-02  -4.40395916e-03  -1.34227730e-03   1.50316724e-02   1.14689633e-01   9.74802430e-02   8.16953368e-02   6.42347352e-02   1.14863690e-01]
[ -2.48934037e-01  -8.34448790e-01   1.20541848e-01  -4.61428969e-02   1.10340611e-02  -1.38274027e-01   1.54718718e-01   1.58408510e-02   1.58655384e-02   4.16141469e-03   3.60515363e-03  -2.61841776e-01  -2.20909661e-01  -1.85873715e-01  -1.43558666e-01   7.68506567e-02]
[ -1.81207641e-01   1.20537511e-01  -9.12297754e-01  -3.37076829e-02   7.93129006e-03  -1.00006758e-01   1.12614657e-01   1.14552887e-02   1.21451662e-02   3.14253621e-03   2.49890970e-03  -1.90783327e-01  -1.60793846e-01  -1.35366758e-01  -1.04288378e-01   5.58544984e-02]
[  7.43758735e-02  -4.61861906e-02  -3.37397718e-02  -9.61551455e-01   1.74836317e-02  -8.98819614e-02  -4.44710999e-02   1.08682106e-02  -1.25474295e-01  -2.41955988e-02   2.48250432e-02   1.14611542e-01   6.31907668e-02   6.82729945e-02  -2.36503319e-04  -5.08569560e-03]
[ -1.25531195e-02   1.09967538e-02   7.90498685e-03   1.74965174e-02  -9.82789579e-01  -1.12302141e-01   9.21095895e-03   1.33248212e-02  -9.61305432e-02  -1.82172436e-02   2.09786819e-02   1.60028081e-02  -1.33975818e-02   8.53570200e-04  -4.18649149e-02   1.82271837e-02]
[  6.38057525e-02   1.24718862e-01   9.00934514e-02   1.04032856e-01   1.19729834e-01  -9.31499844e-01   1.09339291e-01   9.47648922e-02   2.64131010e-02   4.72762859e-03   1.42689857e-01   5.79300210e-02   5.73323636e-02   4.44262653e-02   4.76842316e-02   1.46438318e-01]
[ -2.32816373e-01   1.54719429e-01   1.12619209e-01  -4.44310509e-02   9.24567851e-03  -1.22537590e-01  -8.55441086e-01   1.40088873e-02   2.10810529e-02   5.07876003e-03   2.03310798e-03  -2.46769430e-01  -2.06461090e-01  -1.74497661e-01  -1.32033702e-01   7.09507849e-02]
[ -2.08115816e-02   1.58125336e-02   1.14355020e-02   1.08800825e-02   1.33265812e-02  -8.98683871e-02   1.39824856e-02  -9.89390841e-01  -7.07377866e-02  -1.33519272e-02   1.57628810e-02  -1.76591484e-04  -2.01471347e-02  -7.93638089e-03  -3.77446716e-02   1.70947859e-02]
[ -4.39615929e-03  -1.97439812e-02  -1.49717925e-02   1.27547846e-01   9.66702553e-02   2.64168741e-02  -2.47554289e-02   7.09595290e-02  -9.73213228e-01   5.12198340e-03   1.22666446e-01  -1.26261876e-02  -3.28268619e-03  -6.08434806e-03   6.95146187e-03   6.87239971e-02]
[ -1.34092041e-03  -4.83000977e-03  -3.62959570e-03   2.45018201e-02   1.82691521e-02   4.72827724e-03  -5.70953876e-03   1.33595657e-02   5.12198880e-03  -9.99045890e-01   2.32500551e-02  -2.92806032e-03  -1.07265164e-03  -1.53261064e-03   1.02586508e-03   1.25360928e-02]
[ -3.55257001e-04   3.55897148e-03   2.46608728e-03   2.48380511e-02   2.09761965e-02  -1.32517696e-01   1.99019610e-03   1.57588137e-02  -1.21746730e-01  -2.31431712e-02  -9.73903901e-01   3.63175408e-02  -3.15147533e-03   1.25820724e-02  -4.37034836e-02   1.81394397e-02]
[  1.14687468e-01   2.47573480e-01   1.80328301e-01  -9.55781671e-02  -4.86251696e-03   5.79213957e-02   2.32663340e-01   7.80107025e-03  -1.26341107e-02  -2.92943113e-03  -2.14388161e-02  -8.77792555e-01   1.01651871e-01   8.61957889e-02   6.42372713e-02   9.79062923e-02]
[  9.74807096e-02   2.07531845e-01   1.50995626e-01  -4.62218968e-02   2.31387972e-02   5.73256867e-02   1.93280542e-01   2.67724031e-02  -3.28989498e-03  -1.07390893e-03   1.62177826e-02   1.01654343e-01  -9.13451926e-01   7.24778004e-02   5.72604694e-02   1.03703534e-01]
[  8.16951801e-02   1.75222981e-01   1.27564138e-01  -5.44075733e-02   7.18739633e-03   4.44202822e-02   1.63985606e-01   1.34234104e-02  -6.09057287e-03  -1.53369672e-03  -1.82081989e-03   8.61973527e-02   7.24772812e-02  -9.38929859e-01   4.67971943e-02   7.77429043e-02]
[  6.42376700e-02   1.33212905e-01   9.67155701e-02   1.23871809e-02   4.86178123e-02   4.76808343e-02   1.21890605e-01   4.22879094e-02   6.94643031e-03   1.02498470e-03   5.28282150e-02   6.42416314e-02   5.72628819e-02   4.67996353e-02  -9.59038190e-01   9.34053301e-02]
[ -1.12307485e-01   7.68178288e-02   5.58331499e-02  -5.05776188e-03   1.82415250e-02  -1.46103311e-01   7.09199777e-02   1.71041618e-02  -6.98869555e-02  -1.27692448e-02   1.81593462e-02  -9.49159570e-02  -1.01458672e-01  -7.57289231e-02  -9.22890030e-02  -9.53958411e-01]]
Spectral norm of X =  1.0000572062
Round(X.T * X) =
[[ 1.  0.  0. -0. -0. -0.  0. -0. -0. -0. -0. -0. -0. -0. -0.  0.]
[ 0.  1. -0.  0.  0.  0. -0.  0.  0.  0.  0.  0.  0.  0.  0. -0.]
[ 0. -0.  1.  0.  0.  0. -0.  0.  0.  0.  0.  0.  0.  0.  0. -0.]
[-0.  0.  0.  1. -0. -0.  0. -0. -0. -0. -0. -0. -0. -0. -0.  0.]
[-0.  0.  0. -0.  1. -0.  0. -0. -0. -0. -0. -0. -0. -0. -0.  0.]
[-0.  0.  0. -0. -0.  1.  0. -0. -0. -0. -0. -0. -0. -0. -0.  0.]
[ 0. -0. -0.  0.  0.  0.  1.  0.  0.  0.  0.  0.  0.  0.  0. -0.]
[-0.  0.  0. -0. -0. -0.  0.  1. -0. -0. -0. -0. -0. -0. -0.  0.]
[-0.  0.  0. -0. -0. -0.  0. -0.  1. -0. -0. -0. -0. -0. -0.  0.]
[-0.  0.  0. -0. -0. -0.  0. -0. -0.  1. -0. -0. -0. -0. -0.  0.]
[-0.  0.  0. -0. -0. -0.  0. -0. -0. -0.  1. -0. -0. -0. -0.  0.]
[-0.  0.  0. -0. -0. -0.  0. -0. -0. -0. -0.  1. -0. -0. -0.  0.]
[-0.  0.  0. -0. -0. -0.  0. -0. -0. -0. -0. -0.  1. -0. -0.  0.]
[-0.  0.  0. -0. -0. -0.  0. -0. -0. -0. -0. -0. -0.  1. -0.  0.]
[-0.  0.  0. -0. -0. -0.  0. -0. -0. -0. -0. -0. -0. -0.  1.  0.]
[ 0. -0. -0.  0.  0.  0. -0.  0.  0.  0.  0.  0.  0.  0.  0.  1.]]

I ran the Python script several times until I obtained a matrix $$\rm X$$ that is approximately orthogonal. In other words, the algorithm does not produce such nice results for all choices of $$\rm (a,b,D)$$. Also, the larger the value of $$n$$, the better the results seem to be — which is why I picked $$n=16$$, not $$n=3$$.

Your belief that there should always be a solution is wrong for $n=1$ if $D=a=1$ and $b=-1$, which satisfies all your requirements.

Thus it seems unlikely that one can say anything in general.

Posing the linear and the orthogonality constraint as a least squares problems in the Frobenius norm, and submitting it to an unconstrained optimization routine is probably the only thing one can do.

• What I meant was I believe the input data should always be consistent. I was thinking there is a variation of Gram-Schmidt that can make this work, but I can't seem to figure it out. Nov 19 '12 at 21:15