2
$\begingroup$

The following set of operators $P_1,P_2,P_3,P_4$ and $Q_1,Q_2,Q_3,Q_4$ are related as $P_i=\sum_j w_{ij} Q_j$. How can I find the matrix $W=[w_{ij}]$?

 Q1 = {{0, 0}, {Sqrt[x*y], 0}}
Q2 = {{0, Sqrt[x*(1 - y)]}, {0, 0}}
Q3 = Sqrt[y]*{{Sqrt[1 - x], 0}, {0, 1}}
Q4 = Sqrt[1 - y]*{{1, 0}, {0, Sqrt[1 - x]}}
P1 = {{0, 0}, {Sqrt[x*y], 0}}
P2 = {{0, Sqrt[x*(1 - y)]}, {0, 0}}
P3 = (Sqrt[2 - x - Sqrt[4 - 4*x + x^2 - 4*x^2*y + 4*x^2*y^2]]/
        (Sqrt[2]*Sqrt[1 + (1/4)*
                 ((-x + 2*x*y + 
              Sqrt[4 - 4*x + x^2 - 4*x^2*y + 4*x^2*y^2])/
                      Sqrt[1 - x])^2]))*
     {{-((-x + 2*x*y + Sqrt[4 - 4*x + x^2 - 4*x^2*y + 4*x^2*y^2])/
              (2*Sqrt[1 - x])), 0}, {0, 1}}
P4 =  
     (Sqrt[2 - x + Sqrt[4 - 4*x + x^2 - 4*x^2*y + 4*x^2*y^2]]/
          (Sqrt[2]*Sqrt[1 + (1/4)*
                   ((-x + 2*x*y - 
               Sqrt[4 - 4*x + x^2 - 4*x^2*y + 4*x^2*y^2])/
                        Sqrt[1 - x])^2]))*
       {{-((-x + 2*x*y - Sqrt[4 - 4*x + x^2 - 4*x^2*y + 4*x^2*y^2])/
                (2*Sqrt[1 - x])), 0}, {0, 1}}
$\endgroup$
3
  • $\begingroup$ Do you mean Pi=sum over j for wij*Qj instead of sum over i? Does that mean P1==w11*Q1+w12*Q2+w13*Q3+w14*Q4. Similarly for the other 15 variables? Can you pick a couple of those equations and use Reduce to see if it can find w11 and w12 for example? If that works then maybe you could begin replacing some of the wij with known values. For 16 equations and unknowns and matricies that complicated I'm worried. $\endgroup$
    – Bill
    Sep 25, 2021 at 17:38
  • 1
    $\begingroup$ Thanks, @Bil, just corrected the typo. Actually, as we can see $P_1=Q_1$ and $P_2=Q_2$, so that would mean $w_{12}=w_{13}=w_{14}=w_{21}=w_{23}=w_{24}=0$. $\endgroup$
    – Zubin
    Sep 25, 2021 at 18:07
  • 2
    $\begingroup$ VERY good! Eliminating that many unknowns makes me less worried. It would be nice if there was some high level abstract way of writing your system, but can you try using Reduce on some or all of the manually written equations and see if it can find some more of your unknowns? Maybe after you have a solution you can see a more compact way of writing it. Or maybe someone else can suggest a bright way of doing this. $\endgroup$
    – Bill
    Sep 25, 2021 at 18:32

1 Answer 1

2
$\begingroup$

A quick solution

You can do this in general with standard linear algebra techniques, since the space of matrices forms a vector space under ordinary matrix addition and scalar multiplication! This is made easier to accomplish in Mathematica if we actually just convert our matrices into Mathematica vectors via Flatten.

Then, we could use the following function to extract the coefficients, with vs as the vector(s) we want to represent, and bs the "basis" vectors we want to represent them in terms of:

W[vs_, bs_] := Enclose@ConfirmBy[vs.PseudoInverse[bs], Simplify[#.bs == vs] &]

vs can be a single vector or a list of vectors:

(* Examples: *)
W[{2, 2, -1}, {{1, 1, 0}, {0, 0, 1}}]

(* Out: {2, -1} *)

W[{{2, 2, -1},{0, 0, 5}}, {{1, 1, 0}, {0, 0, 1}}]

(* Out: {{2, -1}, {0, 5}} *)

We need to use PseudoInverse in general in order to handle cases where we don't span the full space. (When we do span the full space, it coincides with Inverse.) We also need Simplify for the symbolic case, otherwise == might wind up unreduced (and therefore not technically True).

So, converting each of our matrices into vectors with Flatten:

w = W[Flatten /@ {P1, P2, P3, P4}, Flatten /@ {Q1, Q2, Q3, Q4}]

This gives us our matrix $w_{ij}$, which can probably be simplified symbolically with some Assumptions on y and x, e.g. Simplify[w, Assumptions -> 1 > x > 0 && 1 > y > 0], or whatever your conditions are!

But we can check that it's at least correct in general, without any assumptions, and calculate $P_i = \sum_{ij} w_{ij}Q_j$:

{P1, P2, P3, P4} == w.{Q1, Q2, Q3, Q4} // Simplify

(* Out: True *)



In general

Setup

Consider any vector $v$ which we wish to express in terms of the vectors $b_1, b_2,\dots$. Suppose we have each $b_i$ expressed in terms of a standard basis, $b_i=\sum_j b_{ij}\hat{e}_j$. We want the coefficients $a_i$ needed to express $v$ as $\sum_i a_i b_i$. But notice that, expressing $b_i$ and $v$ in the standard basis, we equivalently want $\sum v_j \hat{e}_j=\sum_{ij} a_i b_{ij} \hat{e}_j$ and therefore $v_j=\sum_i a_i b_{ij}$ for all $j$. Acting on the right by the inverse of the matrix $b=[b_{ij}]$, denoted $b^{-1}$, we have $$\sum_j v_j (b^{-1})_{jk}=\sum_{ij}a_ib_{ij}(b^{-1})_{jk}=\sum_i a_i\delta_{ik}=a_k$$ In general, though, it might be the case that $\{b_i\}_i$ does not constitute a basis. Either it doesn't span the full space, or it contains linear dependence.

This is nicely solved by PseudoInverse: it exists (and works) even if the vectors are linearly dependent, and can be computed for arbitrary rectangular matrices, both numeric and symbolic. We only need to make sure that whatever we get does actually work; the fact that the PseudoInverse exists for all rectangular matrices necessarily means that sometimes it can't function as a kind of inverse (as sometimes $v$ cannot be expressed in terms of the $b_i$, but the PseudoInverse alone can't keep track of that information for all potential $v$).

Let's make our function a bit more complete in terms of its functionality and its error handling.

Design

Here's how we want our function to work:

  • Allow options (in particular, Assumptions) to be provided which will be passed on to Simplify.

  • Make sure the first argument is either a vector or a list of vectors, and make sure the second is a list of vectors, otherwise send a message and return unevaluated.

  • Compute vs.PseudoInverse[bs], Simplify it according to any Assumptions given, and then check that the result satisfies Simplify[#.bs == vs, Assumptions -> optionalAssumptions] &.

    • If we get False and the vector(s) are definitely not in the span of our "basis" vectors, return a Failure.

    • If the vector might be in the span of our basis matrices depending on symbolic conditions, return a ConditionalExpression.

    • If a solution is found and there's linear dependence among the bs, still provide the solution, but send a message alerting the user that there may be more and, if possible, what they are. MatrixRank and NullSpace seem to struggle with assumptions, e.g. showing linear dependence among symbolic vectors such as {x, y} and {x, s} under the assumption y == s, so we'll use a special function.

(We could also use LinearSolve for this instead of PseudoInverse, and it might be faster; I haven't tried it yet. But PseudoInverse is conceptually nice.)

Code

ClearAll[VectorCoefficients];

Options[VectorCoefficients] = Options[Simplify];

VectorCoefficients::vectmat = 
  "Argument `1` at position `2` is not a matrix or a vector.";

VectorCoefficients::matrix = 
  "Argument `1` at position `2` is not a list of equal-length \
vectors.";

VectorCoefficients::vbsh = 
  "Vector(s) `1` must be the same length as vectors in `2`.";

VectorCoefficients::nlindep = 
  "There is linear dependence among the list of vectors given in \
position `1`; the returned vector(s) of coefficents are unique only \
up to the span of `2`.";

VectorCoefficients::snlindep = 
  "There is linear dependence among the list of vectors given in \
position `1`; the returned vector(s) of coefficents are nonunique at \
least up to the span of the vectors `2`.";

VectorCoefficients::slindep = 
  "There is linear dependence among the list of vectors given in \
position `1`; the solution given may not be unique.";

VectorCoefficients::slindepc = 
  "If `2`, there is linear dependence among the list of vectors given \
in position `1`; in that case (and only that case), the solution \
given may not be unique.";

VectorCoefficients::vfail = 
  "The vector `vector` is not contained in the span of `bvectorlist`.";

VectorCoefficients::vsfail = 
  "Not all vectors in `vectorlist` are contained in the span of \
`bvectorlist`.";

VectorCoefficients[vs_, bs_, opts : OptionsPattern[]] := 
 Module[{sopts = Sequence @@ FilterRules[{opts}, Options@Simplify],
         cs, r, c, lde},
   (* Get the (candidate) coefficients *)

   cs = Simplify[vs . PseudoInverse[bs], sopts];

   (* Get the condition under which the coefficients are valid: *)
   
   r = Simplify[cs . bs == vs, sopts];
   
   (* If r is not literally false, 
      check for linear dependence in the bs and send a Message if there
      is any. This block doesn't do or evalaute to anything else. *)
   
   If[r =!= False,
    (* Produce different messages depending on whether the matrices 
       are both numeric or not. *)
    
    If[AllTrue[bs, NumericQ, 2]
        && (! VectorQ[vs] || AllTrue[vs, NumericQ])
        && (! MatrixQ[vs] || AllTrue[vs, NumericQ, 2]),
     If[MatrixRank[bs] < Length[bs],
      Message[VectorCoefficients::nlindep, 2, NullSpace@Transpose@bs]],
     (* Get the conditions under which symbolic linear dependence
     occurs; if this certainly happens, send a simple message, 
     but if whether it happens depends on the symbols, 
     send a message including the condition. *)
     
     If[MatrixRank[bs] < Length[bs],
      Message[VectorCoefficients::snlindep, 2, NullSpace@Transpose@bs],
      lde = 
       Simplify[
        Reduce[
         Exists[c, c \[Element] Vectors[Length@bs], 
          c != 0 && c . bs == 0], c], sopts];
      If[lde,
       Message[VectorCoefficients::slindep, 2];,
       Null,
       Message[VectorCoefficients::slindepc, 2, lde]]]]];
   
   (* Form the appropriate return value depending on whether r is True, False,
      or nonboolean. *)
   If[r,
    cs,
    If[VectorQ[vs], 
     Failure[
      "SpanMembershipFailure", <|
       "MessageTemplate" :> VectorCoefficients::vfail, 
       "MessageParameters" -> <|"vector" -> vs, "bvectorlist" -> bs|>|>],
     Failure[
      "SpanSubsetFailure", <|
       "MessageTemplate" :> VectorCoefficients::vsfail, 
       "MessageParameters" -> <|"vectorlist" -> vs, "bvectorlist" -> bs|>|>]],
    ConditionalExpression[cs, r]]] /;
  (* Only apply this definition in the first place if the arguments 
     are appropriate; otherwise, leave it unevaluated and send a message *)
  (VectorQ[vs] || MatrixQ[vs] || Message[VectorCoefficients::vectmat, vs, 1])&&
  (MatrixQ[bs] || Message[VectorCoefficients::matrix, bs, 2]) && 
  (Last@Dimensions[vs] == Last[Dimensions[bs]] ||
   Message[VectorCoefficients::vbsh, vs, bs])

Tests/Demos

Here are a bunch of demonstration evaluations you can evaluate in a block as-is. The Echo @ Unevaluated @ ... will show the input expression above each resulting output cell. This should produce some clean output, but also a lot of error messages and Failures, as it's meant to demonstrate that part of the function as well. All messages should be from VectorCoefficients, not from other functions.

Echo@Unevaluated@VectorCoefficients[{{1}, 2}, {{1, 1}, {0, 1}}]

Echo@Unevaluated@VectorCoefficients[{1, 2}, {{{1}, 1}, {0, 1}}]

Echo@Unevaluated@VectorCoefficients[{1, 2}, {{1, 1, 2}, {0, 1, 1}}]

Echo@Unevaluated@VectorCoefficients[{1, 2}, {{1, 1}, {0, 1}}]

Echo@Unevaluated@
  VectorCoefficients[{y, 2, 1}, {{1, 0, 0}, {0, 1, 0}, {x, 0, 1}}]

Echo@Unevaluated@VectorCoefficients[{1, 2}, {{1, 1}, {0, 1}, {0, 1}}]

Echo@Unevaluated@
  VectorCoefficients[{{1, 2}, {3, 4}, {5, 6}}, {{1, 1}, {0, x}}]

Echo@Unevaluated@VectorCoefficients[{1, 2}, {{1, x}, {1, 1}}]

Echo@Unevaluated@VectorCoefficients[{1, 2}, {{0, x}, {0, 1}}]

Echo@Unevaluated@VectorCoefficients[{{1, 2}, {0, 1}}, {{0, x}, {0, 1}}]

Echo@Unevaluated@VectorCoefficients[{2, 2}, {{x, x}, {1, 1}}]

Echo@Unevaluated@
  VectorCoefficients[{2, 2}, {{x, x}, {1, 1}}, 
   Assumptions -> x \[Element] Reals]

Echo@Unevaluated@VectorCoefficients[{y, z}, {{x, 0}}]

Echo@Unevaluated@
  VectorCoefficients[{y, z}, {{x, 0}}, Assumptions -> x != 0]

Echo@Unevaluated@
  VectorCoefficients[{y, z}, {{x, 0}}, 
   Assumptions -> x != 0 && z == 0]

Possible Improvements

  • Allow different return types—for example, Reduce-like conditions telling us what coefficient vectors/matrices can work, or the approximation we get from PseudoInverse instead of the usual Failure—governed by an option.

  • Refactor the testing and handling of linear dependence into a separate function.

  • Provide more complete information on the other solutions due to linear dependence.

  • Allow the information on other solutions due to linear dependence to be provided in the return value, depending on an option.

  • Allow the user to use FullSimplify instead of Simplify, depending on an option (avoided by default because FullSimplify can take longer, but some things may not be caught by Simplify)

  • Allow the user to easily turn off the linear dependence warning via an option

  • When vs is a list of several vectors, consider each one individually, instead of returning a Failure if just one of them is outside the span of the bs.

  • Thread the comparison of the "check" vector(s) and vs and apply And so that we can get better conditions in certain symbolic cases, e.g. ones which aren't {y, 0} == {y, z} but simply z == 0. (For some reason this is unaffected by either Simplify or FullSimplify.

  • Check if LinearSolve is faster or more useful somehow, and understand the differences.

  • Currently VectorCoefficients only handles Mathematica bona-fide vectors and matrices (VectorQ and MatrixQ), but we could potentially expand this to allow linear combinations of arbitrary tensors. Perhaps we would just need to flatten at the appropriate depth (found by comparing the dimensions of the two arguments) when encountering those arrays.

  • We could also generalize the above to allow us to consider the vectors to be at different depths—e.g. instead of a list of basis vectors, use a matrix of basis vectors.

Anyway, this was fun! I was expecting there to be a built-in for handling this task cleanly, honestly, but couldn't find it. Even if there is and I've just somehow missed it, it was still fun! Hope it's helpful! :)

$\endgroup$

Your Answer

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

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