Let's be lazy (and avoid identifier clashes)

```
sage: V=var("v", n=9)
sage: P=matrix(V,nrows=3)
sage: X=matrix(SR,3,3,[[0,1,0],[0,0,1],[1,0,0]])
sage: Pdagger=P.transpose()
sage: Xdagger=X.transpose()
sage: Q=X*Pdagger
sage: B=Q.solve_left(X)
```

Let's be lazy again : don't build a list of equations representing B==X, but a list representing the elements of B-X

```
sage: E=[B[r][c]-X[r][c] for c in range(3) for r in range(3)]
```

Look (lazily) for solutions

```
sage: S=solve(E,V)
```

A loook at the solutions show that they are rational expressions. Are they acceptable ? This entails that their denominators must not be null. But :

```
sage: [[e.denominator().subs(s) for s in S] for e in E]
[[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1]]
```

Check the only acceptable one :

```
sage: [e.subs(S[16]) for e in E]
[0, 0, 0, 0, 0, 0, 0, 0, 0]
```

And, indeed :

```
sage: matrix([v.subs(S[16]) for v in V],nrows=3)==X
True
```

Note that Sympy does not seem to be able to solve this system :

```
sage: sympy.solve(map(lambda e:e._sympy_(), E),map(lambda v:v._sympy_(), V))
[]
```

On the other hand, Mathematica is able to solve the system and does not gives "spurious" solutions. But expressing the exuations is harder in Mathematica than in Sage : that's why I keep the first part in Sage :

```
sage: V=var("v",n=9)
sage: P=matrix(SR,V,nrows=3);P
[v0 v1 v2]
[v3 v4 v5]
[v6 v7 v8]
sage: X=matrix(SR,[[0,1,0],[0,0,1],[1,0,0]]);X
[0 1 0]
[0 0 1]
[1 0 0]
sage: Q=X*(P.transpose())
sage: B=Q.solve_left(X)
sage: E=[B[l][c]==X[l][c] for c in range(3) for l in range(3)];E
[(v2*v7/v1 - v8)*(v0/v1 + v2*(v3 - v0*v4/v1)/(v1*(v2*v4/v1 - v5)))/((v2*v4/v1 - v5)*(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) + v2/(v1*(v2*v4/v1 - v5)) == 0,
-(v0/v1 + v2*(v3 - v0*v4/v1)/(v1*(v2*v4/v1 - v5)))/(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1) == 0,
-(v0/v1 + v2*(v3 - v0*v4/v1)/(v1*(v2*v4/v1 - v5)))*(v4*(v2*v7/v1 - v8)/(v1*(v2*v4/v1 - v5)) - v7/v1)/(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1) + 1/v1 - v2*v4/(v1^2*(v2*v4/v1 - v5)) == 1,
-1/(v2*v4/v1 - v5) - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/((v2*v4/v1 - v5)^2*(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) == 1,
(v3 - v0*v4/v1)/((v2*v4/v1 - v5)*(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) == 0,
(v3 - v0*v4/v1)*(v4*(v2*v7/v1 - v8)/(v1*(v2*v4/v1 - v5)) - v7/v1)/((v2*v4/v1 - v5)*(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) + v4/(v1*(v2*v4/v1 - v5)) == 0,
-(v2*v7/v1 - v8)/((v2*v4/v1 - v5)*(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) == 0,
(1/(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1)) == 1,
(v4*(v2*v7/v1 - v8)/(v1*(v2*v4/v1 - v5)) - v7/v1)/(v6 - (v3 - v0*v4/v1)*(v2*v7/v1 - v8)/(v2*v4/v1 - v5) - v0*v7/v1) == 0]
sage: mathematica.Solve(mathematica.And(*E), V)
{{v0 -> 0, v1 -> 1, v2 -> 0, v3 -> 0, v4 -> 0, v5 -> 1, v6 -> 1, v7 -> 0,
v8 -> 0}}
```

This seems to trigger a bug. I'll need some time to follow it (not so soon).