2
$\begingroup$

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)}}
Improve this question
$\endgroup$
13
  • 1
    $\begingroup$ 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. $\endgroup$
    – Community Bot
    Jun 20, 2022 at 12:09
  • 1
    $\begingroup$ Doesn't A . PseudoInverse[M] suit your needs? $\endgroup$
    – J. M.'s lack of A.I. ♦
    Jun 20, 2022 at 14:07
  • 1
    $\begingroup$ You might wish to include an example where LeastSquares[]/PseudoInverse[] is unable to give a solution with nonnegative entries (if you have one). $\endgroup$
    – J. M.'s lack of A.I. ♦
    Jun 20, 2022 at 14:27
  • 1
    $\begingroup$ 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]. $\endgroup$
    – J. M.'s lack of A.I. ♦
    Jun 20, 2022 at 15:01
  • 1
    $\begingroup$ @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. $\endgroup$
    – Farid Shahandeh
    Jun 20, 2022 at 15:11

2 Answers 2

Reset to default
3
$\begingroup$

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 )

Improve this answer
$\endgroup$
4
  • $\begingroup$ 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.) $\endgroup$
    – J. M.'s lack of A.I. ♦
    Jun 20, 2022 at 14:19
  • $\begingroup$ 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. $\endgroup$
    – Farid Shahandeh
    Jun 20, 2022 at 14:25
  • $\begingroup$ @FaridShahandeh No I don't think so. Both variants only give solution with minimal Norm[X] $\endgroup$
    – Ulrich Neumann
    Jun 20, 2022 at 14:27
  • $\begingroup$ @UlrichNeumann Yes. The norm does not guarantee nonnegativity. I now remember why about a year ago I used FindInstance in my heuristic code. $\endgroup$
    – Farid Shahandeh
    Jun 20, 2022 at 14:44
2
$\begingroup$

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?

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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