We use the following pseudo Newton iteration, that gives very quickly an orthogonal matrix.
$V_{k+1}=0.5(V_k+V_k^{-T})$.
Indeed, let $f:V\rightarrow V^TV-I_n$ and $Df_V:H\rightarrow H^TV+V^TH$. The standard Newton iteration is $V_{k+1}=V_k-(Df_{V_k})^{-1}(f(V_k))$, but $Df_{V_k}$ is not invertible; yet $H=0.5(V_k-V_k^{-T})$ is a solution of $Df_{V_k}(H)=f(V_k)$. For this solution $V_{k+1}=0.5(V_k+V_k^{-T})$.
A trial is composed with the following $3$ steps.
Step 1. We choose an element of $[L,U]$: $V_0=kL+(1-k)U$ where $k$ has the form $i/10,i=0,\cdots,10$.
Step 2. We use the above iterative algorithm with initial point $V_0$ and with $10$ iterations.
Step 3. If the obtained limit is in the box, then we stop. Otherwise, we continue with the next value of $k$.
The time of calculation for $1000=10^3$ trials is $17"$. There are on average $5$ failures for $10000=10^4$ trials.
Below, a copy of the procedure in Maple (the result is in $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;
EDIT. Case 1. There is some solutions in the box and we want to find an approximation of one such solution. A Maple specialist advised me to use the command NPSolve; unfortunately, NPSolve only works (in general) if you give it an initial point not too far from a solution. Yet we can do as follows:
STEP 1. As in my above procedure, we pick some point $V0=kL+(1−k)U$ where k has the form $i/10,i=0,⋯,10$ on the segment $[L,U]$ and we project it in $V_1$ on $O(3)$. In a few cases, none of the $11$ points $V_1$ is in the box; yet, some are very close to the edge of the box.
STEP 2. We use NLPSolve (this soft accepts only the closed boxes) with $V_1$ as initial point; then the software uses the tangent vector space of $O(3)$ in $V_1$ and, with a gradient method, pushes this point towards the edge of the box. In general we obtain by this method a point $V_2$ that is very close to the edge of the box and we are done. For 20000=2.10^4
trials (as in my above procedure), I did not suffer any failure (on the $11$ points $V_1$, NLPSolve always works for many of them)!!
Of course there is no proof that the method is without fail. One must even be able to find a counter-example by making sure that there is only one solution (on the edge of the box). Of course, the OP did not take a lot of risk by offering a big bonus for rigorous proof. He will be forgiven because 500000 dollars is the price of a nice house.
CASE 2. There are no solutions in the box and we want to prove this fact. Except if you find a condition in the form $x[i,j]>1$, I think we cannot avoid a formal proof. In these conditions, I do not see any other methods than the Gröbner's basis method. Over $\mathbb{C}$, there are black boxes that do the job when the number of relations is not too large; unfortunately, it is known that in general the problem is much more complicated in the real case; in particular, it is necessary to use gradient methods. Anyway, I don't know how to do the job with Maple.
