Given two full-rank matrices $A \in \mathbb{R}^{n \times p}, B \in \mathbb{R}^{n \times k}$ and vectors $y \in \mathbb{R}^n, u \in \mathbb{R}^p, v \in \mathbb{R}^k$ I'd like to solve an optimization problem of the form

$$ \arg\min_{u, v} \|y - (A u) \circ (Bv)\|^2$$

where $\circ$ denotes the Hadamard product (aka elementwise product) and the norm is the euclidean vector norm.

I can find the solution by alternate minimization (iteratively fixing one variable and solving for the other one), but it is terribly slow.In practice I have a lot of these problems where the $y$ varies but $A, B$ are always the same so ideally I would like to express the solution in terms of a matrix factorization of $A$ and/or $B$.

My question is: has this problem been addressed somewhere in the literature? Is it possible to find a closed form solution in terms of a matrix factorization of $A, B$ analogous to the (ordinary) least squares?

  • $\begingroup$ The dude's name is HadamarD. $\endgroup$ – Mariano Suárez-Álvarez Feb 13 '14 at 9:17
  • $\begingroup$ Have you tried standard solvers or packages? Have you tried Gauss-Newton for nonlinear least squares? I haven't tried to compute the jacobian of $(Au)*(Bv)$, but it shouldn't be that complicate, as the function is bilinear. $\endgroup$ – user251257 Jul 2 '15 at 1:31
  • $\begingroup$ Note that what I am asking for are either pointers to the literature where this problem has been addressed or a closed form solution. $\endgroup$ – Fabian Pedregosa Jul 6 '15 at 15:33
  • 1
    $\begingroup$ I suspect there might be no closed-form solution for this problem since the problem isn't convex and there doesn't seem to be a nice transformation (at least at a first glance) which would make it so. Perhaps your best bet might be to just to alternating optimization, as you mentioned. $\endgroup$ – Guillermo Angeris Aug 22 '17 at 14:40

Here I'd change $\alpha \in \mathbb{R}^p, \ \beta \in \mathbb{R}^k$ just for notation coherence.

One way to transform this problem into an ordinary least squares (OLS) is:

We have our main problem with Hadamard product:

$$ \begin{equation} \text{argmin}_{\alpha,\beta}\|y - \bar{y} \|_2^2, \ \ \ (1) \end{equation} $$

were $\bar{y} = (A\alpha \odot B\beta)$ is our approximation vector.

We can rewrite each row of our approximation vectos as:

$$ \bar{y}_i = \left( A_i \alpha \right)\left(B_i \beta\right), $$

where $A_i$ and $B_i$ denotes the $i^\text{th}$ row of the corfesponding matrix.

Just changing the notation, we have $$ \begin{align} \bar{y}_i & = \left(\begin{array}{l l l} a_{i,1} B_i, \ldots, a_{i,p} B_i \end{array} \right) \left( \begin{array}{} \alpha_1 \beta \\ \vdots \\ \alpha_p \beta \end{array} \right), \end{align} $$

$$ (A\alpha \odot B\beta)_i = (A_i \otimes B_i)(\alpha \otimes \beta), $$

where $\otimes$ denotes the kronecker product.

Using this equality we can define our new design matrix as $ H_i = (A_i \otimes B_i )$, and solution vector as $x = (\alpha \otimes \beta)$, then we can solve the OLS:

$$ \text{argmin}_{x}\|y - Hx\|_2^2, $$

finding our global solution, due the convexity in $x$.

In order to retrieve the solution vectors, can rewrite $x$ as a matix, for instance, product of our solution vectors $\alpha$ and $\beta$,

$$ X = \left( \begin{array}{lcr} \alpha_1 \beta_1 & \ldots & \alpha_1 \beta_k \\ \alpha_2 \beta_1 & \ldots & \alpha_2 \beta_k\\ \vdots & \ddots & \vdots \\ \alpha_p \beta_1 & \ldots & \alpha_p \beta_k \end{array} \right) = \alpha \beta^T $$

We use the singular value decomposotion (SVD) to represent our matrix, as

$$ X=U\Sigma V^T, $$

as $Rank(X)=1$, then using the first left and right singular vectors $u$ and $v$, and the first singular value $\sigma$,

$$ \bar{\alpha} = \gamma\sqrt\sigma u, \ \ \text{and} \ \ \bar{\beta} = \frac{1}{\gamma}\sqrt\sigma v, $$

where $\gamma$ depends of the scale of $\alpha$ and $\beta$, and also the bias, $|\gamma|>0$.

The vectors solutions are scale sensitive, so the solution is not unique and additional assumptions should be made.

Hope this will be usefull

  • $\begingroup$ Thanks for your answer, but you are assuming that the least squares solution on $H$ followed by the best rank-1 approximation to that solution will be the global solution, which is generally not the case. More formally, you are assuming that $\text{argmin}\|y - X (\alpha \otimes \beta)\|$ = $\text{argmin}\|X^{\dagger} y - (\alpha \otimes \beta)\|$ which is in general not true (take for example H = [[1, 1, 1, 1], [1, 0, 0, 0]], y = [1, 2], with your solution the residual is 0.23 but the solution (1, 1), (-2, 1.5) has zero residual) $\endgroup$ – Fabian Pedregosa Mar 4 '14 at 8:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.