Maximize $l_1$ norm with unitary matrix Given an invertible matrix $A \in \mathbb{C}^{n\times n}$. How to find
$$
U^* = \max_{\text{$U$ with $U^H U = I$}} \lVert U A\rVert_1,
$$
where $\lVert\cdot\rVert_1$ is the entrywise 1-norm, i.e., $\lVert A\rVert_1 = \sum_{i,j} \lvert A_{ij}\rvert$ and $\cdot^H$ denotes the complex conjugate transpose?
For $A = I$, the solution is any complex Hadamard matrix, e.g., a scaled discrete Fourier matrix. 
 A: Too long to comment:
You might want to try out the following sub-optimal approaches, as the original problem is likely to be in NP-class (citations needed!). I am also going to assume that $A$ has only real entries, for simplicity.
1) Take a look at the paper "Orthogonalization of vectors with minimal adjustment" in Biometrika. The first optimality criterion might be a good surrogate for your cost function. Intuitively, it finds the an orthogonal basis which has least deviation (in terms of dot product) from the a given set of vectors (here columns of $A$).
2) You might also want to investigate the sub-optimality of $Q$ where $A=QR$ (QR decomposition). Changing the order in which Gram-Schmidt is done will yield different QR decompositions, and so you might also want to use this degree of freedom judiciously.
3) Lastly, the greedy way. Intuitively, suppose one were to find a unit vector $q_1$ such that $||q_1A||_1$ is maximized. Post this, one can find a unit vector maximizing $||q_2A||_1$, subject to the additional constraint $q_1q^\top_2=0$. This can be done till one obtains a orthonormal matrix. What remains now is:
$$
\max_{q\in R^{n\times 1}} ||qA||_1~\mbox{subject to}~qq^\top = 1 ~\&~ qC = 0.
$$
Intuitively, a unit vector has highest 1-norm if all its entries are equal in magnitude, or its alignment along on of the vectors $[\pm 1,\cdots,  \pm 1]$ is maximum. This can be used in the following iterative routine. Choose any unit vector $q^{(0)}$ such that $q^{(0)}C=0$. Run the following convex routine iteratively till convergence (there is no guarantee of convergence):
$$
q^{(k+1)} = \arg \max_{qA\in R^{n\times 1}} q\left({\mbox{sign}(q^{(k)})}^\top\right) ~\mbox{subject to}~ qq^\top\leq 1 ~\&~ qC=0.
$$
(4) One could also try the following simpler idea. Note that 
$\begin{bmatrix} I & X\\X^\top &I \end{bmatrix} \succeq 0 $ is a convex relaxation to orthonormal constraint (using Schur complement). And a the maximum of a linear cost functional with this relaxed constraint will always yeild an optima at its boundary, which is the set of orthonormal matrix. Using the same logic as in point (3), one can write the code below:
import numpy as np
import cvxpy as cvx
X = cvx.Variable((N,N))
Q,R = np.linalg.qr(np.random.randn(N,N))
Q = np.sign(Q@A)
count = 0
while(count<=20):
constraints = [cvx.vstack((cvx.hstack((np.eye(N),X)),cvx.hstack((X.T,np.eye(N))))) >> 0]

prob = cvx.Problem(cvx.Maximize(cvx.trace((X @ A)@(Q.T))), constraints)

prob.solve()

Q = np.sign((X.value)@A)

count = count + 1

print(np.sum(np.abs(X.value @ A)))
Z,R = np.linalg.qr(np.random.randn(N,N))
print(np.sum(np.abs(Z@A)))
Hope this helps.
A: Here an alternative solution, which outperformed @DSM solution (4) in all tested cases.
Without loss of generality, we assume that $\lVert A \rVert_F = 1$. The optimal value for the $\ell_1$ norm is attained by the unitary Hadamard matrix $H$, e.g., DFT matrix.
As a proxy cost function, we use therefore
$$
\hat{U} = \min_U \| |U A| - |H| \|_F, 
$$
where $|\cdot|$ is the element-wise absolute values. This is equivalent for an optimal set of phases $\hat{P}$ with $|\hat{P}_{ij}| = 1$ such that 
$$
\hat{U} = \min_U \| U A - |H| \circ \hat{P} \|_F,
$$
where $\circ$ denotes the element-wise (Hadamard) product. This can be solved iteratively by unitary Procrustes solution such that
$$
U^{(i+1)} = \min_U \| U A - |H| \circ P^{(i)} \|_F \\
P^{(i+1)} = U^{(i+1)}A \oslash |U^{(i+1)}A|, 
$$
where $\oslash$ is the element-wise (Hadamard) division. These iterations are guaranteed to converge: 
Because of the Procrustes solution is the global minimum for the Frobenius norm, we have
$$
\| U^{(i+1)} A - |H| \circ P^{(i)} \|_F \leq \| U^{(i)} A - |H| \circ P^{(i)} \|_F
$$
Then, updating the phase also reduces the error
$$
\| U^{(i+1)} A - |H| \circ P^{(i+1)} \|_F \leq \| U^{(i+1)} A - |H| \circ P^{(i)} \|_F,
$$
which is essentially a element-wise version of
$$
\phi = \min_\theta \left(a e^{\imath \phi} - b e^{\imath \theta} \right)^2,
$$
where $a$, $b$, $\phi$, $\theta$ are real valued.
