Cannot solve Ax=b for A using python if there is no unique solution

Question

I am trying to solve a linear system of equations of the form Ax=b where A is a square matrix which is unknown, while the vectors x and b are known. I am interested in a physics problem where my A is a unitary matrix, which rotates a random unit vector x to b which is of the form [1,0,0...0]. One cannot find a unique solution to A. I just need one of the many possible solutions. Since I want to scale it to larger dimensions, I want the computer to solve it.

I tried following the answer to the question in this link , but by A was not unitary.

I am thinking of using cvxpy to give me a solution, and I wrote a code taking the dimension to be three for a start. I don't know if cvxpy is suitable when there is no unique solution. The code is below. It gives an error in the constraint part saying "0-dimensional array given. Array must be at least two-dimensional". I think what I have is a valid representation of a matrix in cvxpy. I am happy to be corrected.

Can someone help me with this problem? If you have another way of doing it, that is also fine.

import numpy as np
import cvxpy as cp
n = 3
alpha=0.01
np.random.seed(1)
b= [1,0,0]
x = np.random.randn(n)+1j*np.random.randn(n)
x=x/np.linalg.norm(x)


#Construct the problem.
A = cp.Variable((n,n),complex=True)
objective = cp.Minimize(cp.sum_squares(A @ x - b))
constraints = [A.H==np.linalg.inv(A)]
prob = cp.Problem(objective, constraints)

# The optimal objective value is returned by `prob.solve()`.
result = prob.solve()
print(A.value)

Traceback:

Traceback (most recent call last):
    File "/Users/sreerampg/untitled1.py", line 22, in <module>     constraints = [A.H==np.linalg.inv(A)]
    File "<__array_function__ internals>", line 6, in inv    File "/Applications/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py", line 539, in inv     _assert_stacked_2d(a)
    File "/Applications/anaconda3/lib/python3.7/site-packages/numpy/linalg/linalg.py", line 197, in _assert_stacked_2d     
    'at least two-dimensional' % a.ndim)  LinAlgError: 0-dimensional array given. Array must be at least two-dimensional

Answer 1

It seems like you could approach this analytically, since you're looking for a rotation. In 3 dimensions, you could directly construct the rotation matrix by computing the angle between x and b, and the vector perpendicular to both:

import numpy as np
from scipy.spatial.transform import Rotation

def find_rot(a, b):
    a /= np.linalg.norm(a)
    b /= np.linalg.norm(b)
    angle = np.arccos(np.dot(a, b))
    vec = np.cross(a, b)
    return Rotation.from_rotvec(vec * angle / np.linalg.norm(vec))

x = np.random.rand(3)
b = np.array([1., 0., 0.])
A = find_rot(x, b)
print(x, b, A.apply(x), A.as_matrix(), sep='\n')

In higher dimensions it's not as direct since you can't just use the cross product, although there's probably a way to generalize the above using exterior product for which I'd need to brush up on my linear algebra skills more.

Barring that, you can construct a series of "higher dimensional rotation matrices" that successively zero out each dimension, of the form:

[1 . . . . 
 . 1 
 .   1
 .     . 
 .       cosT -sinT
 .       sinT cosT]

Where T is the angle theta between your vector x and the unit vector [0, 0, ...., 1, 0], then repeat moving down one dim at a time. Essentially you "rotate out" your vector from each dimension, then multiply them all together to get your final unitary matrix.

Answer 2

You want Ax = b, where A is unitary, so this is equivalent to A^* b = x, where A^* is the complex transpose (Hermitian). Look for A^* instead. Since b = [1,0,...], the vector b is the first unit vector, so you are looking for unitary matrices with x on the first column (and this determines the first row):

A^* = [x1, x2^*,...,xn^*|
       x2, ...

in other word, the solution is a matrix A with x^* (the complex conjugate of x) on the first column, the rest of the components of x on the first row, and a unitary matrix of size (n-1)x(n-1) filling the rest of the matrix.

(As you see, for this to have a solution, x1, the first component of x must be real)