# solve -- solve linear equation(s)

## Synopsis

• Usage:
X = solve(A,B)
• Inputs:
• A, , or , 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 => , default value false, whether to use the least squares method, in the case when the ring is RR or CC
• MaximalRank => , 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.000449264 seconds i31 : time X = solve(A,B, MaximalRank=>true); -- used 0.000094387 seconds i32 : norm(A*X-B) o32 = 5.11185069084045e-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.118574 seconds i39 : time X = solve(A,B, MaximalRank=>true); -- used 0.112653 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 = | -5.551e-16 | | -3.331e-16 | | 0 | 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.