Here is the matlab code I used:
function mathoverlow_triangular
global n k nk
n=10;
k=4;
nk=n*k;
options=optimset('Jacobian','on');%'Display','iter',,'DerivativeCheck','on'
exitflag=0;
while exitflag~=1 %until a solution has been found
[x,~,exitflag]=fsolve(@(x) f(x),rand(2*nk,1),options);
end
[~,~,A]=f(x)
%find 0-1 pattern
nz=double(abs(A)>=abs(A'));
save('sol','nz');
for reps=1:1e4 %try this many times finding an integer solution
p=randperm(n);
ind=p(1:k);
B1=nz;
pm=1-2*(rand(n,k)>.5);
is=1;
B1(:,ind)=nz(:,ind).*pm;%.*randi(4,n,k)
for j=k+1:n
nu=null(B1(B1(:,p(j))==0,ind),'r');
%exclude those which give 0 on diagonal
if size(nu,2)>0
nu=nu(:,B1(p(j),ind)*nu~=0);
end
if size(nu,2)==0
is=0;
break;
else
B1(:,p(j))=B1(:,ind)nu(:,randi(size(nu,2)));
end
end
if is
A=B1;
[N,D]=rat(A);
A=round(A.(ones(n,1)*lcm_array(D)));
A=round(A./(ones(n,1)*gcd_array(A)));
A=round(A./(gcd_array(A')'*ones(1,n)));
if max(abs(A(:)))==1
%write tex code
for i=1:n
stri='';
for j=1:n-1
stri=[stri num2str(A(i,j)) ' & '];
end
disp([stri num2str(A(i,end)) '\']);
end
disp('');
save('sol','A');
end
end
end
end
function [b]=lcm_array(A)
if size(A,1)==1
b=A;
else
b=lcm(lcm_array(A(1:end-1,:)),A(end,:));
end
end
function [b]=gcd_array(A)
if size(A,1)==1
b=A;
else
b=gcd(gcd_array(A(1:end-1,:)),A(end,:));
end
end
function [res,J,UV] = f(x)
global n k nk
U=reshape(x(1:nk),[n k]);
V=reshape(x(nk+1:end),[n k]);
UV=UV';
Tnk=transposeT(n,k);
res=UV'.UV-eye(n);
Tnn=transposeT(n,n);
n2=n^2;
dUV=sparse(1:n2,1:n2,UV(:),n2,n2);
UVt=UV';
dUVt=sparse(1:n2,1:n2,UVt(:),n2,n2);
J=(dUVt+dUVTnn)[kron(V,speye(n)) kron(speye(n),U)*Tnk];
end
function T = transposeT(n,k)
%derivative of transpose map
nk=n*k;
u=reshape(1:nk,[n k]);
v=u';
T=sparse(u(:),v(:),ones(nk,1),nk,nk);
end
function mathoverlow_triangular
global n k nk
n=10;
k=4;
nk=n*k;
options=optimset('Jacobian','on');%'Display','iter',,'DerivativeCheck','on'
exitflag=0;
while exitflag~=1 %until a solution has been found
[x,~,exitflag]=fsolve(@(x) f(x),rand(2*nk,1),options);
end
[~,~,A]=f(x)
%find 0-1 pattern
nz=double(abs(A)>=abs(A'));
save('sol','nz');
for reps=1:1e4 %try this many times finding an integer solution
p=randperm(n);
ind=p(1:k);
B1=nz;
pm=1-2*(rand(n,k)>.5);
is=1;
B1(:,ind)=nz(:,ind).*pm;%.*randi(4,n,k)
for j=k+1:n
nu=null(B1(B1(:,p(j))==0,ind),'r');
%exclude those which give 0 on diagonal
if size(nu,2)>0
nu=nu(:,B1(p(j),ind)*nu~=0);
end
if size(nu,2)==0
is=0;
break;
else
B1(:,p(j))=B1(:,ind)*nu(:,randi(size(nu,2)));
end
end
if is
A=B1;
[N,D]=rat(A);
A=round(A.*(ones(n,1)*lcm_array(D)));
A=round(A./(ones(n,1)*gcd_array(A)));
A=round(A./(gcd_array(A')'*ones(1,n)));
if max(abs(A(:)))==1
%write tex code
for i=1:n
stri='';
for j=1:n-1
stri=[stri num2str(A(i,j)) ' & '];
end
disp([stri num2str(A(i,end)) '\\']);
end
disp('');
save('sol','A');
end
end
end
end
function [b]=lcm_array(A)
if size(A,1)==1
b=A;
else
b=lcm(lcm_array(A(1:end-1,:)),A(end,:));
end
end
function [b]=gcd_array(A)
if size(A,1)==1
b=A;
else
b=gcd(gcd_array(A(1:end-1,:)),A(end,:));
end
end
function [res,J,UV] = f(x)
global n k nk
U=reshape(x(1:nk),[n k]);
V=reshape(x(nk+1:end),[n k]);
UV=U*V';
Tnk=transposeT(n,k);
res=UV'.*UV-eye(n);
Tnn=transposeT(n,n);
n2=n^2;
dUV=sparse(1:n2,1:n2,UV(:),n2,n2);
UVt=UV';
dUVt=sparse(1:n2,1:n2,UVt(:),n2,n2);
J=(dUVt+dUV*Tnn)*[kron(V,speye(n)) kron(speye(n),U)*Tnk];
end
function T = transposeT(n,k)
%derivative of transpose map
nk=n*k;
u=reshape(1:nk,[n k]);
v=u';
T=sparse(u(:),v(:),ones(nk,1),nk,nk);
end