Solution to a linear equation with minimum rank

Question

Suppose I have a linear equation, $X.M=A$, for given matrices $M$ and $A$. I can use FindInstance to produce a solution $X$.

How can I add the constraint that $X$ has the least rank?

Indeed, I could add Det[$X$]==0 in FindInstance when $X$ is square, and iterate. But, $X$ may be non-square.

Example 1:

Suppose

M={{0, 1, 0, 0, 1, 1}, {0, 1, 0, 1, 1, 0}, {0, 1, 1, 0, 0, 1}, {0, 1, 1,
   1, 0, 0}, {1, 0, 0, 0, 1, 1}, {1, 0, 0, 1, 1, 0}, {1, 0, 1, 0, 0, 
  1}, {1, 0, 1, 1, 0, 0}}

and

A={{1, 0, 1/2, 1/2, 1/2, 1/2}, {0, 1, 1/2, 1/2, 1/2, 1/2}, {1/2, 1/2, 1,
   1/2, 0, 1/2}, {1/2, 1/2, 1/2, 1, 1/2, 0}, {1/2, 1/2, 0, 1/2, 1, 1/
  2}, {1/2, 1/2, 1/2, 0, 1/2, 1}}

Then,

X1 = Array[f, {Dimensions[[A]], Dimensions[M][[1]]}]; 
X1 /. FindInstance[
   X1 . M == A, 
   Flatten[X], NonNegativeReals][[1]]

returns

X1={{0, 0, 0, 0, 0, 1/2, 1/2, 0}, {0, 1/2, 1/2, 0, 0, 0, 0, 0}, {0, 0, 1/
  2, 0, 0, 0, 0, 1/2}, {0, 1/2, 0, 0, 0, 0, 0, 1/2}, {1/2, 0, 0, 0, 0,
   1/2, 0, 0}, {0, 0, 1/2, 0, 1/2, 0, 0, 0}};
X//MatrixRank
6

There is, however, a solution of rank four, namely,

X2={{0, 0, 0, 0, 1/2, 0, 0, 1/2}, {0, 1/2, 1/2, 0, 0, 0, 0, 0}, {0, 0, 1/
  2, 0, 0, 0, 0, 1/2}, {0, 1/2, 0, 0, 0, 0, 0, 1/2}, {0, 1/2, 0, 0, 1/
  2, 0, 0, 0}, {0, 0, 1/2, 0, 1/2, 0, 0, 0}};
X2//MatrixRank
4

How can I enforce FindInstance to look for a solution with the lowest rank, i.e., returns X2 rather than X1?

The only thing I got to work is

FindInstance[
       X1 . M == A&&Det[X1]==0, 
       Flatten[X], NonNegativeReals]

But, this works only for square matrices. In general, I encounter non-square ones.

As suggested in the comments, LeastSquares provides a way to find the lowest rank solution. However, as in the example above, it's important that the matrix entries of the solution are nonnegative. Here is an example showing the difficulty with the LeastSquates.

Example 2:

Let

M={{0, 0, 1, 1}, {0, 1, 1, 0}, {1, 0, 0, 1}, {1, 1, 0, 0}}

and

A={{1, 0, 0, 1}, {1, 1, 0, 0}, {0, 1, 1, 0}, {0, 0, 1, 1}}

Then, LeastSquares returns solutions with negative entries:

X={{1/4, -(1/4), 3/4, 1/4}, {-(1/4), 1/4, 1/4, 3/4}, {1/4, 3/4, -(1/4), 
  1/4}, {3/4, 1/4, 1/4, -(1/4)}}

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking. — Community, Jun 20 '22 at 12:09
Hi Farid, welcome to Mathematica.SE, please take the [tour] so you familiarize yourself with the expectations on this site. Here it's considered helpful and polite to show your own efforts and share your data and code attempts in a well-formatted form form, so we can quickly see the problem you are facing. It is NOT clear at all that this is a Wolfram Mathematica programming question. Probably you are looking for math.stackexchange.com? Please help us help you, and edit your question to offer a minimum working example of code. — rhermans, Jun 20 '22 at 12:29
@rhermans Hi, and thanks for your comment. I've used the mathematica.SE for a long time, it's just the first time I'm posting a question. So, I assume I'm familiar with the basics. Unfortunately, on this occasion, I'm not able to share my code, as it's too big and my question concerns only a few lines of it. I'm not sure using FindInstance or constraining a linear equation has anything to do with math.SE. I believe it's straightforward to choose two matrices M and A and see if there is a way to find the solution X with minimum rank (when the solution is not unique). — Farid Shahandeh, Jun 20 '22 at 12:37
@UlrichNeumann I thought about that. But, I couldn't find a way to implement it, because X is a priori unknown. — Farid Shahandeh, Jun 20 '22 at 13:54
@J.M. Thanks both. I'm trying to fit this into my actual code and check. It's important for me that the entries of the solution are all nonnegative too (and thus the NonNegativeReals in my example). — Farid Shahandeh, Jun 20 '22 at 14:22
You might wish to include an example where LeastSquares[]/PseudoInverse[] is unable to give a solution with nonnegative entries (if you have one). — J. M.'s missing motivation, Jun 20 '22 at 14:27
@J.M. I added an example to show the case with negativities. — Farid Shahandeh, Jun 20 '22 at 14:42
In the particular case you gave, it seems like a solution with nonnegative entries is not possible: m = {{0, 0, 1, 1}, {0, 1, 1, 0}, {1, 0, 0, 1}, {1, 1, 0, 0}}; a = {{1, 0, 0, 1}, {1, 1, 0, 0}, {0, 1, 1, 0}, {0, 0, 1, 1}}; pinv = PseudoInverse[Transpose[m]]; ns = First[NullSpace[Transpose[m]]]; Reduce[Thread[Flatten[pinv . Transpose[a] + KroneckerProduct[ns, ConstantArray[t, Length[a]]]] > 0], t]. — J. M.'s missing motivation, Jun 20 '22 at 15:01
@J.M. There is one nonnegative solution, X1={{0, 0, 1, 0}, {0, 0, 0, 1}, {0, 1, 0, 0}, {1, 0, 0, 0}} which is of rank 4. — Farid Shahandeh, Jun 20 '22 at 15:11
I found the solution, with the help of Ulrich Neuman and J.M.'s slightly less busy's comments above, to be the nonnegative least squares (NNLS). The algorithm is not built-in to Mathematica. But, there is one available online. I've posted it separately here. — Farid Shahandeh, Jun 21 '22 at 09:17

Ulrich Neumann · Answer 1 · 2022-06-20T14:14:24.850

3

Your transposed equation looks like Transpose[M].Transpose[X]==Transpose[A]

LeastSquares finds the minimal solution

X = Transpose[LeastSquares[Transpose[M], Transpose[A]]]
(*{{0, 0, 0, 0, 1/4, 1/4, 1/4, 1/4}, 
  {1/4, 1/4, 1/4, 1/4, 0, 0, 0,0},
  {0, 0, 1/4, 1/4, 0, 0, 1/4, 1/4},
  {0, 1/4, 0, 1/4, 0, 1/4, 0, 1/4},
  {1/4, 1/4, 0, 0, 1/4, 1/4, 0, 0}, 
  {1/4, 0, 1/4, 0, 1/4, 0, 1/4,0}}*)

Hope it helps!

This solution this coincides with A.PseudoInverse[M] (thanks @ J. M.'s slightly less busy comment )

edited Jun 20 '22 at 14:14

answered Jun 20 '22 at 14:08

Ulrich Neumann

53,729
2
23
55

In the exact case, the snippet I gave might make the intent clearer. If the given matrices had inexact entries, however, your solution is certainly preferable. (This is similar to why one might prefer to use LinearSolve[] over Inverse[] in evaluating a certain matrix expression.) – J. M.'s missing motivation Jun 20 '22 at 14:19
Yes, it might be the case that the inputs have inexact values, unlike my example. Is there a way to guarantee the nonnegativity of the entries of the solution, given the nonnegativity of M and A? I'm checking if that's the case for the least-square problem in general. – Farid Shahandeh Jun 20 '22 at 14:25
@FaridShahandeh No I don't think so. Both variants only give solution with minimal Norm[X] – Ulrich Neumann Jun 20 '22 at 14:27
@UlrichNeumann Yes. The norm does not guarantee nonnegativity. I now remember why about a year ago I used FindInstance in my heuristic code. – Farid Shahandeh Jun 20 '22 at 14:44

Alex Trounev · Answer 2 · 2022-06-26T02:52:35.467

We can use Nonnegative Least Squares Algorithm (NNLS) discussed here to solve matrix equation X.M=A with constraints Table[X[[i,j]]>=0,{i,m},{j,n}]. First we define functions

bitsToIndices[v_] := 
 Module[{l = Select[Table[i, {i, Length[v]}], v[[#]] == 1 &]}, l]; 
compZ[aa_, bb_, ff_] := 
 Module[{Ap, Q, R}, 
  Ap = Transpose[Transpose[aa][[bb]]]; {Q, R} = QRDecomposition[Ap]; 
  Inverse[R] . Q . ff];
NNLS1[A_, f_] := 
  Module[{x, zeroed, w, t, z, q, \[Alpha], i, zeroedSet, positiveSet, 
    toBeZeroed}, zeroedSet := bitsToIndices[zeroed];
   positiveSet := bitsToIndices[1 - zeroed];
   x = 0 A[[1]];
   zeroed = 1 - x;
   w = Transpose[A] . (f - A . x);
   Do[If[zeroedSet != {} && Max[w[[zeroedSet]]] > 0,
      t = Position[w zeroed, Max[w zeroed], 1, 1][[1]][[1]];
      zeroed[[t]] = 0;
      z = 0 x;
      z[[positiveSet]] = compZ[A, positiveSet, f];
      Do[If[Min[z] < 0,
         \[Alpha] = Infinity;
         Do[
          If[zeroed[[q]] == 0 && 
             z[[q]] < 0, \[Alpha] = 
              Min[\[Alpha], x[[q]]/(x[[q]] - z[[q]])];
            ];, {q, 1, Length[x]}];
         x = x + \[Alpha] (z - x);
         toBeZeroed = Select[positiveSet, Abs[x[[#]]] < 10^-13 &];
         zeroed[[toBeZeroed]] = 1;
         x[[toBeZeroed]] = 0;
         z = 0 x;
         z[[positiveSet]] = compZ[A, positiveSet, f];, 
         Break[]];, {Infinity}
       ]; x = z;
      w = Transpose[A] . (f - A . x);, Break[]];, {Infinity}
    ]; Return[x];];

Then we transform matrix equation to the vector form as follows

X = Array[u,Dimensions[Transpose[M]]]; var = 
 Flatten[X, 2];eq = Flatten[X . M - A, 1]; {v, mat} = CoefficientArrays[Table[eq[[i]] == 0, {i, Length[eq]}], var];

Finally we solve equation mat.var==-v with using NNLS1 and transform solution into matrix form

s = NNLS1[mat // Normal, -v // Normal];Xsol = X /. Table[var[[i]] -> s[[i]], {i, Length[var]}];

In a case of Example 2 we have

M = {{0, 0, 1, 1}, {0, 1, 1, 0}, {1, 0, 0, 1}, {1, 1, 0, 0}};
A = {{1, 0, 0, 1}, {1, 1, 0, 0}, {0, 1, 1, 0}, {0, 0, 1, 1}};

Solution s={0, 0, 1, 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 0} and consequently

Xsol={{0, 0, 1, 0}, {0, 0, 0, 1}, {0, 1, 0, 0}, {1, 0, 0, 0}}

Test that X.M-A=0 and

Xsol // MatrixRank
Out[]= 4

It could be useful to test this code with rectangular matrix as well. But in a case of Example 1 we have

M = {{0, 1, 0, 0, 1, 1}, {0, 1, 0, 1, 1, 0}, {0, 1, 1, 0, 0, 1}, {0, 
    1, 1, 1, 0, 0}, {1, 0, 0, 0, 1, 1}, {1, 0, 0, 1, 1, 0}, {1, 0, 1, 
    0, 0, 1}, {1, 0, 1, 1, 0, 0}};
A = {{1, 0, 1/2, 1/2, 1/2, 1/2}, {0, 1, 1/2, 1/2, 1/2, 1/2}, {1/2, 
    1/2, 1, 1/2, 0, 1/2}, {1/2, 1/2, 1/2, 1, 1/2, 0}, {1/2, 1/2, 0, 
    1/2, 1, 1/2}, {1/2, 1/2, 1/2, 0, 1/2, 1}}; 
Xsol={{0, 0, 0, 0, 1/2, 0, 0, 1/2}, {1/2, 0, 0, 1/2, 0, 0, 0, 0}, {0, 0, 1/
  2, 0, 0, 0, 0, 1/2}, {0, 1/2, 0, 0, 0, 0, 0, 1/2}, {1/2, 0, 0, 0, 0,
   1/2, 0, 0}, {1/2, 0, 0, 0, 0, 0, 1/2, 0}};

This is exactly what we got with using FindInstance since Xsol // MatrixRank is 6 in this case. Therefore, we try Code 2 from my answer here

NNLSFindMinimum[A_, f_] := 
  Module[{nbx = Length[First[A]], xi, x, axf, xinit}, 
   xi = Array[x, nbx];
   axf = A . xi^2 - f;
   xinit = PseudoInverse[A] . f;
   If[And @@ (# >= 0 & /@ xinit), xinit, 
    fm = FindMinimum[Evaluate[axf . axf], 
      Evaluate[Sequence @@ Transpose[{xi, xinit}]], 
      MaxIterations -> 1000];
    xi^2 /. fm[[2]]]];
s2 = NNLSFindMinimum[m // Normal, -v // Normal]
(Out[]= {0, 0, 0, 0, 1/4, 1/4, 1/4, 1/4, 1/4, 1/4, 1/4, 1/4, 0, 0, 

0, 0, 0, 0, 1/4, 1/4, 0, 0, 1/4, 1/4, 0, 1/4, 0, 1/4, 0, 1/4, 0, 1/4, 

1/4, 1/4, 0, 0, 1/4, 1/4, 0, 0, 1/4, 0, 1/4, 0, 1/4, 0, 1/4, 0})
X2 = X /. Table[var[[i]] -> s2[[i]], {i, Length[var]}]
(Out[]= {{0, 0, 0, 0, 1/4, 1/4, 1/4, 1/4}, {1/4, 1/4, 1/4, 1/4, 0, 0,
   0, 0}, {0, 0, 1/4, 1/4, 0, 0, 1/4, 1/4}, {0, 1/4, 0, 1/4, 0, 1/4, 
  0, 1/4}, {1/4, 1/4, 0, 0, 1/4, 1/4, 0, 0}, {1/4, 0, 1/4, 0, 1/4, 0, 
  1/4, 0}})
X2 . M - A
(Out[]= {{0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 
  0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}, {0, 0, 0, 0, 0, 0}})
X2 // MatrixRank
(Out[]= 4)

But how can we prove that Code 2 by Jean - Claude Poujade minimize rank as well?

Solution to a linear equation with minimum rank

2 Answers2

Linked