Python precision and performance on large integer matrix products

Question

I want to calculate the product of two quantities A and B and then compute the product modulo F. In general, A and B are matrices with large numbers so the operation is the conventional matrix multiplication. Let's consider for simplicity A and B as scalars. The Python code with all the methods I have tested is the following:

from __future__ import division
import numpy as np
from sympy import Matrix
from sympy import *

A = 2251875000001
B = 28839630

F = 33232924804801

#Method 1: pure Python
C1 = (A*B) % F

A_mat = np.matrix(A).astype(np.int64)
B_mat = np.matrix(B).astype(np.int64)

#Method 2: numpy matrices and pure Python
C2 = (A_mat*B_mat) % F

#Method 3: numpy
C3 = np.dot(A, B) % F

#Method 4: numpy
C4 = np.dot(A_mat, B_mat) % F

#Method 5: Python objects and numpy
A_mat2 = np.matrix(A, dtype=object)
B_mat2 = np.matrix(B, dtype=object)
C5 = np.dot(A_mat2, B_mat2) % F
C5 = np.concatenate(C5).astype(np.int64)

#Method 6: sympy
A_sp = Matrix(1,1,[A])
B_sp = Matrix(1,1,[B])
C6 = A_sp*B_sp
f = lambda x: x%F
C6 = C6.applyfunc(f)
print(C6)
C6 = matrix2numpy(C6, dtype='int64')

Now the result should be

(2251875000001 * 28839630) mod 33232924804801 = 25112458407047

When I run the above code what I get for C1 == 25112458407047 is correct for this example, but when I test multiplication of large matrices most of the entries I am getting with this method are wrong. However, the values of C2, C3, C4 all equal to 12062945446480 which is wrong. C5 and C6 are also calculated correctly.

You can assume that 64bit precision of int64 is more than enough for the numbers I am working with, but the default 32bit is not.

I tried Sympy (see method 6) since it's supposed to have arbitrary precision according to here.

I am using Python 2.7.14, Numpy 1.14.2 and Sympy 1.1.1.

My first question is why do I get some of the results wrong in the above code?

Finally, methods 5 and 6 even though they are always correct, they seem to be slower than the other methods. Can you explain that? How do the above methods differ in terms of complexity of matrix multiplication and what do you suggest?

EDIT

I thought this should be clear from the functions I use in the code, but anyway, the operation I am interested in is conventional matrix multiplication. Also, I apologize for the fact that my example was actually an overflow, but the question of performance is still valid.

`A*B <= 2**64 -1` returns `False` for me, meaning 64 bit precision is actually not enough to store the result of multiplication of A and B. — eozd, May 18 '18 at 15:22
You could use this [rule](https://math.stackexchange.com/questions/2416119/rules-for-modulus-and-multiplication) to avoid `A` and `B` multiplication — Brenlla, May 18 '18 at 15:25
It is confusing when you say you want to calculate the product of two quantities A and B, but then you say the operation is the conventional matrix multiplication. Matrix multiplication also involves [sums](https://en.wikipedia.org/wiki/Matrix_multiplication#Definition). If you want the element-wise product, you can use `np.arrays` as follows: `C=((A%F)*(B%F))%F` — Brenlla, May 18 '18 at 15:59
@Brenlla, in this case elementwise multiplication gives the same results as the matrix multiplication. The limitation is the size of `A*B` using `int64`. Using `float` or `float128` arrays gives the right value within 9-12 significant figures. — hpaulj, May 18 '18 at 16:24
@hpaulj Yes, that's true. In the question it is mentioned that A and B are matrices, but in the example A and B are scalars for simplicity. So I understand in the real code A and B are matrices, but unclear if we need to apply matrix or element-wise multiplication — Brenlla, May 18 '18 at 16:30
With `object` dtype arrays, the values are Python integers, and the math is done as in case 1. But it has to 'iterate', so will be slower. — hpaulj, May 18 '18 at 16:33

score 0 · Accepted Answer · 2018-05-18T19:50:09.433

0

As noted in comments, A*B does not fit in 64 bits, so it gets truncated, making the int64-based computation return wrong results. Since int64 is not enough to hold A*B, the next best thing is to use arrays with data type "object" (essentially your method 5, though it's time to leave np.matrix behind):

A_obj = np.array(A).astype(object)
B_obj = np.array(B).astype(object)
np.dot(A_obj, B_obj) % F

returns 25112458407047.

It's slower than int64 multiplication, but still faster than SymPy (by the factor of about 15 in my test).

The existence of np.matrix is a remnant of dark ages of mimicking Matlab syntax; it's not recommended in new code. In my test using np.matrix was about 3 times as slow.

methods 5 and 6 seem to be slower than the other methods.

Handling more complex objects takes longer. It's one thing to process an array of int64 numbers with a low-level routine that is capable of multiplying integers; it's another to process an array of objects that have whatever kind of multiplication operators implemented on them. The former can be (and is) done at C or Fortran level, returning the results to Python; but to handle Python objects, one has to stay in Python all the time.

edited May 18 '18 at 19:50

answered May 18 '18 at 16:39

In method 5, I am using `concatenate` since this data is supposed to be sent from one machine to another via MPI, specifically MPI4py, and the object data type is not supported. So, from what you are telling me I understand that the fastest and reliable way to do matrix product is `A*B` or `np.dot(A,B)` if the numbers are small and if the numbers are large requiring high precision we need to convert to Python objects. – mgus May 18 '18 at 16:58
for `np.matrix`, `dot` and `*` are the same thing. What confused me is the 1st sentence "I want to calculate the product of two quantities A and B" – Brenlla May 18 '18 at 19:12
Good point; I forgot all about np.matrix`. Replacing it by `np.array` resulted in about 3x speedup. @mgus – May 18 '18 at 19:53
@user4412195 Did you notice this speedup even when you use Python `object` as the datatype before the multiplication? Can you tell me what you tried? Because, remember I need this kind of precision since I was thinking but `np.int64` was enough for me but it's not. – mgus May 18 '18 at 20:03
@mgus Yes, I used same datatype for matrix and array. But it does not matter what I see in my notebook; `timeit`in your environment with real data. – May 18 '18 at 22:27

Python precision and performance on large integer matrix products

1 Answers1