solve -- solve linear equation(s)

Synopsis

Usage:

X = solve(A,B)
Inputs:
- A, a matrix, or a mutable matrix, of size m by n over a field R, which can be one of: ZZ/p, GF(p^n), QQ, RR, CC
- B, of the same type of matrix as A, over the same ring, of size m by r
Optional inputs:
- ClosestFit => a Boolean value, default value false, whether to use the least squares method, in the case when the ring is RR or CC
- MaximalRank => a Boolean value, default value false, declares to the system that the matrix is full rank. In some cases, this can dramatically speed up the computation. If the matrix is not full rank, then the results are potentially meaningless.
- Invertible (missing documentation) => ..., default value false,
- Precision => ..., default value 0
Outputs:
- X, of the same type of matrix as A, over the same ring, such that $AX=B$

Description

i1 : kk = ZZ/101;

i2 : A = matrix"1,2,3,4;1,3,6,10;19,7,11,13" ** kk

o2 = | 1  2 3  4  |
     | 1  3 6  10 |
     | 19 7 11 13 |

              3       4
o2 : Matrix kk  <-- kk

i3 : b = matrix"1;1;1" ** kk

o3 = | 1 |
     | 1 |
     | 1 |

              3       1
o3 : Matrix kk  <-- kk

i4 : x = solve(A,b)

o4 = | 2  |
     | -1 |
     | 34 |
     | 0  |

              4       1
o4 : Matrix kk  <-- kk

i5 : A*x-b

o5 = 0

              3       1
o5 : Matrix kk  <-- kk

i6 : kk = GF(25)

o6 = kk

o6 : GaloisField

i7 : a = kk_0

o7 = a

o7 : kk

i8 : A = matrix"a,a+1,a+2,3a,4;a-1,1,2a,6,10;19,7,a,11,13" ** kk

o8 = | a   a+1 a+2 -2a -1 |
     | a-1 1   2a  1   0  |
     | -1  2   a   1   -2 |

              3       5
o8 : Matrix kk  <-- kk

i9 : b = matrix"1;-a+1;1" ** kk

o9 = | 1    |
     | -a+1 |
     | 1    |

              3       1
o9 : Matrix kk  <-- kk

i10 : x = solve(A,b)

o10 = | -a    |
      | -2a+1 |
      | -2a   |
      | 0     |
      | 0     |

               5       1
o10 : Matrix kk  <-- kk

i11 : A*x-b

o11 = 0

               3       1
o11 : Matrix kk  <-- kk

i12 : kk = QQ

o12 = QQ

o12 : Ring

i13 : A = matrix"1,2,3,4;1,3,6,10;19,7,11,13" ** kk

o13 = | 1  2 3  4  |
      | 1  3 6  10 |
      | 19 7 11 13 |

               3       4
o13 : Matrix QQ  <-- QQ

i14 : b = matrix"1;1;1" ** kk

o14 = | 1 |
      | 1 |
      | 1 |

               3       1
o14 : Matrix QQ  <-- QQ

i15 : x = solve(A,b)

o15 = | -7/47  |
      | 54/47  |
      | -18/47 |
      | 0      |

               4       1
o15 : Matrix QQ  <-- QQ

i16 : A*x-b

o16 = 0

               3       1
o16 : Matrix QQ  <-- QQ

Over RR_{53} or CC_{53}, if the matrix A is non-singular and square, then highly optimized lapack routines will be called.

i17 : printingPrecision = 4;

i18 : A = matrix "1,2,3;1,3,6;19,7,11" ** RR

o18 = | 1  2 3  |
      | 1  3 6  |
      | 19 7 11 |

                 3         3
o18 : Matrix RR    <-- RR
               53        53

i19 : b = matrix "1;1;1" ** RR

o19 = | 1 |
      | 1 |
      | 1 |

                 3         1
o19 : Matrix RR    <-- RR
               53        53

i20 : x = solve(A,b)

o20 = | -.1489 |
      | 1.149  |
      | -.383  |

                 3         1
o20 : Matrix RR    <-- RR
               53        53

i21 : A*x-b

o21 = | 2.22e-16  |
      | -2.22e-16 |
      | 0         |

                 3         1
o21 : Matrix RR    <-- RR
               53        53

i22 : norm oo

o22 = 2.22044604925031e-16

o22 : RR (of precision 53)

i23 : clean(1e-15, A*x-b)

o23 = | 0 |
      | 0 |
      | 0 |

                 3         1
o23 : Matrix RR    <-- RR
               53        53

If you know that your matrix is square, and invertible, then providing the hint: MaximalRank=>true allows Macaulay2 to choose the fastest routines. For small matrix sizes, it should not be too noticeable, but for large matrices, the difference in time taken can be dramatic.

i24 : printingPrecision = 4;

i25 : N = 40

o25 = 40

i26 : A = mutableMatrix(CC_53, N, N); fillMatrix A;

i28 : B = mutableMatrix(CC_53, N, 2); fillMatrix B;

i30 : time X = solve(A,B);
     -- used 0.000278346 seconds

i31 : time X = solve(A,B, MaximalRank=>true);
     -- used 0.000046598 seconds

i32 : norm(A*X-B)

o32 = 5.5121203545145e-15

o32 : RR (of precision 53)

Over higher precision RR or CC, these routines will be much slower than the lower precision lapack routines.

i33 : N = 100

o33 = 100

i34 : A = mutableMatrix(CC_100, N, N); fillMatrix A;

i36 : B = mutableMatrix(CC_100, N, 2); fillMatrix B;

i38 : time X = solve(A,B);
     -- used 0.142103 seconds

i39 : time X = solve(A,B, MaximalRank=>true);
     -- used 0.117551 seconds

i40 : norm(A*X-B)

o40 = 1.49157827468970981408235588593e-28

o40 : RR (of precision 100)

Giving the option ClosestFit=>true, in the case when the field is RR or CC, uses a least squares algorithm to find a best fit solution.

i41 : kk = RR_53;

i42 : A = matrix"1,2,3,4;1,3,6,10;19,7,11,13" ** kk

o42 = | 1  2 3  4  |
      | 1  3 6  10 |
      | 19 7 11 13 |

                 3         4
o42 : Matrix RR    <-- RR
               53        53

i43 : b = matrix"1;1;1" ** kk

o43 = | 1 |
      | 1 |
      | 1 |

                 3         1
o43 : Matrix RR    <-- RR
               53        53

i44 : x1 = solve(A,b, ClosestFit=>true)

o44 = | -.1899 |
      | .6399  |
      | .3367  |
      | -.275  |

                 4         1
o44 : Matrix RR    <-- RR
               53        53

i45 : A*x1-b

o45 = | -6.661e-16 |
      | -6.661e-16 |
      | -3.553e-15 |

                 3         1
o45 : Matrix RR    <-- RR
               53        53

Giving both options ClosestFit and MaximalRank allows Macaulay2 to call a faster algorithm.

i46 : x2 = solve(A,b, ClosestFit=>true, MaximalRank=>true)

o46 = | -.1899 |
      | .6399  |
      | .3367  |
      | -.275  |

                 4         1
o46 : Matrix RR    <-- RR
               53        53

i47 : A*x2-b

o47 = | 0         |
      | 0         |
      | 3.109e-15 |

                 3         1
o47 : Matrix RR    <-- RR
               53        53

Caveat

(1) This function is limited in scope, but has been designed to be much faster than generic algorithms. (2) If the matrix is a square invertible matrix, giving the option MaximalRank=>true can strongly speed up the computation. (3) For mutable matrices, this function is only currently implemented for densely encoded matrices.

Ways to use solve :

solve(Matrix,Matrix)
solve(MutableMatrix,MutableMatrix)
solve(Matrix,Vector) (missing documentation)

For the programmer

The object solve is a method function with options.