# How to optimize an approximated matrix multiplication?

[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?

• I'll point out that the condition for invertibility is not correct. The matrix $I-\alpha X$ is invertible if $\alpha<\frac{1}{\lambda_\text{max}(X)}$. Commented Mar 13, 2023 at 22:21
• What does "is selected from identity" mean? Separately, I am confused by the statement that the matrix may have repeated columns but is also a permutation matrix. That can't be right. A permutation matrix is invertible; a matrix with repeated columns is not. Commented Mar 13, 2023 at 22:24
• Thank you for the helpful comments! And my apologies for the misdescriptions. $X$ cannot be a permutation matrix if it has repeated columns. I re-edited my question and present an example when some of the columns are repeated. And yes, the invertibility of $I-\alpha X$ should meet condition $\alpha<1/\lambda_{\max}(X)$, but in this specific case $\alpha<1$ would be fine I think. Commented Mar 14, 2023 at 3:13
• I think I still don't quite understand what "is a standard basis" means. Are you saying that the columns of the matrix are unit coordinate vectors (with one one and the rest being zero)? Commented Mar 14, 2023 at 15:37
• If something is a combinatorial problem, no reformulation is going to change that. If your problem is large (i.e., if the size of the matrices is large), you will not be able to solve the problem exactly. Commented Mar 16, 2023 at 16:02

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.
• Thank you! This makes the question more explicit and easy for me now! However, I still have a small question according to the answer, if I may ask. How to compute matrix $B$? My thought is to maximize the value function by choosing a $X$, the selection of $X$ is by random sampling its columns with probability $p_k$, and the mean gives the approximated maxima. If I do this by random sampling, all I need to do is computing at most $n$ probabilities $p_k$ for $k=1,\cdots,n$ and let the mean of sample does the rest. But under the transformed case, the specifications of $B$ is unknown. Commented Mar 14, 2023 at 3:21
• Yes, the reformulation I have here does not take into account your constraints. I only realized that you had some constraints after writing most of the answer. What one would need to do -- but I'm not sure this is possible here -- is to translate what the constraints on $X$ mean for $B$. Commented Mar 14, 2023 at 15:40