How to optimize an approximated matrix multiplication?

Question

[UPDATING]

The old one is a simplified version of the current one. Here is a solution based on the answer proposed by professor Bangerth down below. To describe what I am trying to do, first rewrite the question into $$\max_X\|(I-\alpha X)^{-1}X\circ AE\|_F^2$$ where $X,A$ still are $n$ by $n$ square matrices, $E$ is a vector with all the elements are one such that $E=(1,\cdots,1)^\top$. According to professor Bangerth's transformation, the objective becomes $$\max_B\left\|\frac{1}{\alpha}(B-I)\circ AE\right\|_F^2$$ where $B=(I-\alpha X)^{-1}$. Now the gradient of objective with respect to $B$ is $$\nabla_B\left(\frac{1}{\alpha^2}\text{Tr}(B-I)\circ(AEE^\top A^\top)\circ(B-I)^\top\right) =\frac{1}{\alpha^2}(B-I)\circ AE(I\circ AE)^\top $$ By the stochastic gradient descent (SGD) algorithm, $B$ follows $$B^{t+1} = B^t + \eta \nabla_B$$ So initialize $B$ as $B^0 = (I-\alpha X^0)^{-1}$ for a guess $X^0$, $B$ can keep iterating till convergence.

The constraint of $B$ is the specification of matrix $X$, which is a sparse matrix and each column has and only has one one. To update $X$ for each iteration, compute the objective, denoted by $J$, with the $B^t$ from SGD, find the location of maxima in each row of $$A + \alpha J$$ then mark the same location in $X$ as one, the rest keeps being zero, this called $X'$. In Python, this can be done by

ini  =  np.argmax(A + alpha * J, axis = 1)
X    =  np.zeros([n,n])
X[np.arange(n),ini]  =  1

The updated $B'$ can be obtained by $X'$ such that $B'=(I-\alpha X')^{-1}$. We can also put $B'$ into SGD to get $B^{t+1}$, these iterations stop when $\|B'-B^t\|<\epsilon$.

The whole thing can be implemented theoretically and pratically, I think. But I still have one pivotal question:

How to add the constraint to the objective or express it as a panelty function?

---

[OLD]

Suppose the objective I try to maximize is $$\max_{X} \|(I - \alpha X)^{-1}XA\|_F$$ where $X$ is the matrix needs to be pinned down, $\alpha$ is a scalar, and $\|\cdot\|_F$ is the Frobenius norm. Note that $(I-\alpha X)^{-1}$ is invertible only when $\alpha\neq1$, so assume $\alpha<1$. All matrices are $n$ by $n$ square matrices. An extra condition of matrix $X$ is that the vectors $X(:,k)$ for $k=1,\cdots,n$ is a standard basis vector selected by the optimization, and $X$ may have repeated columns, for example a 3 by 3 $X$ could be

$$ X = \left( \begin{matrix} 1&0&0\\ 0&1&1\\ 0&0&0 \end{matrix} \right) $$

One possible way to obtain the maxima of the objective is to apply the matrix multiplication via random sampling, where the objective is rewritten as

$$\max_{p_k}\left\| \frac{1}{n}\sum_{k=1}^n\frac{1}{p_k}B(:,k)X(k,:)A\right\|_F$$

where $B$ is the inverse and is approximated by Neumann series such that $$B = (I - \alpha X)^{-1} = \sum_{i=0}^n (\alpha X)^k$$

Note that we do not know the form or the specification of $X$ in advance but choose each column of $X$ by random sampling with probability $p_k$.

Now I have 3 questions:

1. How to compute $B$ when it is a random sampling approximation?

2. Is this method too complex to implement numerically?

3. Is there any better way to solve this puzzle?

Answer 1

This is a linear least squares problem if you just look at it the right way. Write $$ B = (I-aX)^{-1}, $$ then $X = \frac{1}{a}(I-B^{-1})$ and $$ (I-aX)^{-1}XA = B\frac{1}{a}(I-B^{-1})A = \frac{1}{a}(B-I)A. $$ The optimization problem then becomes $$ \max_B \left\|\frac{1}{a}(B-I)A\right\|_F^2, $$ where I have used that the maximizer is the same whether you take the square of the objective function or not. In any case, the objective function is now a square function of $B$, for which the optimality conditions are a linear function of $B$.

The form does suggest that without further constraints on $B$, the objective function grows with $B$ and so no maximizer exists.