# Solution to a Hadamard product least squares

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?

• The dude's name is HadamarD. – Mariano Suárez-Álvarez Feb 13 '14 at 9:17
• 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. – user251257 Jul 2 '15 at 1:31
• Note that what I am asking for are either pointers to the literature where this problem has been addressed or a closed form solution. – Fabian Pedregosa Jul 6 '15 at 15:33
• 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. – 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

• 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) – Fabian Pedregosa Mar 4 '14 at 8:45