Matriz ortogonal con restricciones de caja

Question

Constantes dadas $\ell, u \in \mathbb{R}^{3 \times 3}$ y el siguiente sistema de restricciones en $P \in \mathbb{R}^{3 \times 3}$ $$ P^T P = I_{3 \times 3},\quad \ell_{ij} \leq P_{ij} \leq u_{ij}, $$ Me gustaría encontrar una matriz $P$ que satisfaga este sistema, o determinar que es inviable. ¿Existe una forma computacionalmente eficiente de realizar esta tarea?

La solución no tiene que ser de forma cerrada, pero debe ser un algoritmo implementable en un ordenador que se ejecute rápidamente, explotando el hecho de que es un $9$ problema dimensional.

Un algoritmo numérico que converge a una solución también es una muy buena opción, pero debe haber una prueba de que efectivamente converge a una solución del problema.

Answer 1

Tenemos un sistema de ecuaciones cuadráticas y desigualdades lineales en $\mathrm X \in \mathbb R^{3 \times 3}$

$$\begin{aligned} \rm X^\top X = \rm I_3\\ \rm L \leq X \leq \rm U\end{aligned}$$

donde $\rm L \leq U$ Por supuesto. El casco convexo de la grupo ortogonal $\mathrm O (3)$ se define por $\mathrm X^\top \mathrm X \preceq \mathrm I_3$ o, de forma equivalente, por la desigualdad $\| \mathrm X \|_2 \leq 1$ . Tenemos entonces el siguiente problema de viabilidad (convexo)

$$\begin{array}{ll} \text{minimize} & 0\\ \text{subject to} & \rm L \leq X \leq \rm U\\ & \| \mathrm X \|_2 \leq 1\end{array}$$

Para evitar la solución de la matriz cero, busquemos una solución en el límite de la región factible (convexa). Para ello, generamos una matriz aleatoria $\mathrm C \in \mathbb R^{3 \times 3}$ y minimizar $\rm \langle C, X \rangle$ en cambio

$$\begin{array}{ll} \text{minimize} & \langle \mathrm C, \mathrm X \rangle\\ \text{subject to} & \rm L \leq X \leq \rm U\\ & \| \mathrm X \|_2 \leq 1\end{array}$$

---

Experimento numérico

Supongamos que queremos encontrar un $3 \times 3$ matriz que es ortogonal, cuyas entradas en la diagonal principal están en $[-0.95, 0.95]$ y cuyas entradas fuera de la diagonal principal están en $[-0.5, 0.5]$ .

Uso de NumPy para generar aleatoriamente una matriz $\rm C$ y CVXPY para resolver el programa convexo,

from cvxpy import *
import numpy as np

n = 3

C = np.random.rand(n,n)

I = np.identity(n)
J = np.ones((n,n))

L = -0.95 * I + (-0.5 * (J - I))
U =  0.95 * I + ( 0.5 * (J - I))

X = Variable(n,n)

# define optimization problem
prob = Problem( Minimize(trace(C.T * X)), [ L <= X, X <= U, 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)

que produce lo siguiente

-2.04689293246
optimal
X = 
[[-0.94511744  0.20404426 -0.2650367 ]
 [-0.30732554 -0.84243613  0.44563059]
 [ 0.13860874 -0.499968   -0.85597357]]
Spectral norm of X =  1.00521536292
Round(X.T * X) = 
[[ 1. -0. -0.]
 [-0.  1. -0.]
 [-0. -0.  1.]]

He ejecutado el script unos cuantos (quizás $5$ ) hasta obtener una matriz $\rm X$ que es (aproximadamente) ortogonal. En otras palabras, el algoritmo no produce resultados tan buenos para todo opciones de $\rm C$ .

Answer 2

Utilizamos la siguiente pseudo iteración de Newton, que da muy rápidamente una matriz ortogonal.

$V_{k+1}=0.5(V_k+V_k^{-T})$ .

En efecto, dejemos que $f:V\rightarrow V^TV-I_n$ y $Df_V:H\rightarrow H^TV+V^TH$ . La iteración estándar de Newton es $V_{k+1}=V_k-(Df_{V_k})^{-1}(f(V_k))$ pero $Df_{V_k}$ no es invertible; sin embargo $H=0.5(V_k-V_k^{-T})$ es una solución de $Df_{V_k}(H)=f(V_k)$ . Para esta solución $V_{k+1}=0.5(V_k+V_k^{-T})$ .

Un juicio se compone de lo siguiente $3$ pasos.

Paso 1. Elegimos un elemento de $[L,U]$ : $V_0=kL+(1-k)U$ donde $k$ tiene la forma $i/10,i=0,\cdots,10$ .

Paso 2. Utilizamos el algoritmo iterativo anterior con el punto inicial $V_0$ y con $10$ iteraciones.

Paso 3. Si el límite obtenido está en la casilla, nos detenemos. En caso contrario, continuamos con el siguiente valor de $k$ .

El tiempo de cálculo para $1000=10^3$ los juicios es $17"$ . Hay una media de $5$ fallos por $10000=10^4$ ensayos.

A continuación, una copia del procedimiento en Maple (el resultado está en $V$ ).

restart;
with(LinearAlgebra);
 Digits := 20;

roll := rand(1 .. 3); 

n := 3; id := IdentityMatrix(n):

 ve := Vector(11);

 vv := Vector([.5, .6, .4, .7, .3, .8, .2, .9, .1, 1, 0]);

 t:=time(): 
 for l from 1 to 10000 do

  K:=RandomMatrix(n,n,shape=skewsymmetric):

  U:=evalf((id-K).(id+K)^(-1)): 

 Z1:=Matrix(3,3,(i,j)->0.3*roll()):Z2:=Matrix(3,3,(i,j)->-0.1*roll()): 

  U1:=U+Z1:U2:=U+Z2:

test := [10, 10]; co := 0;
while test<>[0,0] and co<11 do  

co:=co+1:k:=vv[co]:

    V:=k* U1+(1-k)* U2:  
for i to 10 do 

V := .5*(V+1/Transpose(V)) end do;

test:=[0,0]: 

R1:= U+Z1-V:R2:=U+Z2-V: 

 for i from 1 to 3 do 

 for j from 1 to 3 do  

if R1[i,j]<0  then test:=test+[1,0]: fi:  

 if R2[i,j]>0  then test:=test+[0,1]: fi:   

 od:od: 

ve[co]:=test: 

if co=11  and test<>[0,0] then print(ve[7..11],U,Z1,Z2); fi: 

  od:  od:  

time()-t;

EDITAR. Caso 1. Hay algunas soluciones en la caja y queremos encontrar una aproximación de una de esas soluciones. Un especialista de Maple me aconsejó utilizar el comando NPSolve; desgraciadamente, NPSolve sólo funciona (en general) si se le da un punto inicial no muy alejado de una solución. Sin embargo, podemos hacer lo siguiente:

PASO 1. Como en mi procedimiento anterior, elegimos algún punto $V0=kL+(1−k)U$ donde k tiene la forma $i/10,i=0,⋯,10$ en el segmento $[L,U]$ y lo proyectamos en $V_1$ en $O(3)$ . En algunos casos, ninguno de los $11$ puntos $V_1$ está en la caja; sin embargo, algunos están muy cerca del borde de la caja.

PASO 2. Utilizamos NLPSolve (este soft acepta sólo las cajas cerradas) con $V_1$ como punto inicial; entonces el software utiliza el espacio vectorial tangente de $O(3)$ en $V_1$ y, con un método de gradiente, empuja este punto hacia el borde de la caja. En general obtenemos por este método un punto $V_2$ que está muy cerca del borde de la caja y ya está. Para 20000=2,10^4 ensayos (como en mi procedimiento anterior), no sufrí ningún fallo (en el $11$ puntos $V_1$ , ¡¡NLPSolve siempre funciona para muchos de ellos!!)

Por supuesto, no hay pruebas de que el método no tenga fallos. Incluso hay que poder encontrar un contraejemplo asegurándose de que sólo hay una solución (en el borde de la caja). Por supuesto, el OP no se arriesgó mucho al ofrecer una gran bonificación por una prueba rigurosa. Se le perdonará porque 500000 dólares es el precio de una bonita casa.

CASO 2. No hay soluciones en la caja y queremos demostrar este hecho. Excepto si se encuentra una condición de la forma $x[i,j]>1$ Creo que no podemos evitar una prueba formal. En estas condiciones, no veo más métodos que el de la base de Gröbner. En $\mathbb{C}$ En la actualidad, hay cajas negras que hacen el trabajo cuando el número de relaciones no es demasiado grande; por desgracia, se sabe que en general el problema es mucho más complicado en el caso real; en particular, es necesario utilizar métodos de gradiente. De todos modos, no sé cómo hacer el trabajo con Maple.

Answer 3

He aquí una sugerencia. Puede eliminar el $P^TP=I$ de la siguiente manera.

$P^TP=I$ significa que $P$ es una matriz de rotación. Toda matriz de rotación puede escribirse como $\exp([\omega])$ , donde $\exp$ es la matriz exponencial y $$[\omega]=\left( \begin{array}{ccc} 0 & -a_3 & a_2 \\ a_3 & 0 & -a_1 \\ -a_2 & a_1 & 0 \\ \end{array} \right)$$ es una matriz simétrica sesgada. Además, a partir de aquí ( Matriz exponencial de una matriz simétrica sesgada ) se puede encontrar la matriz exponencial analíticamente como $$x = \sqrt{a_1^2+a_2^2+a_3^2}; \quad P=\exp([\omega])= I + \dfrac{\sin x}{x}[\omega] + \dfrac{1-\cos x}{x^2}[\omega]^2.$$

Lo que se me ocurre es realizar una búsqueda por fuerza bruta en $a_1$ , $a_2$ y $a_3$ para ver si cualquier valor para ellos resulta en un $P$ que satisfaga sus restricciones de límites.