Menu
▾
▴
Re: [Maxima-discuss] Error with large matrices
|
From: Oscar B. <osc...@gm...> - 2024-08-11 13:02:21
|
The matrices involved here are sparse and have a simple patterned
structure. It is likely that there is a better approach than using explicit
calculations with large matrices.
If you do want to use explicit calculations with large matrices then you
will definitely want to use something that is designed for sparse matrices
(my impression from previous replies is that Maxima does not have this).
For example if this matrix defines a set of linear equations then you
should use an iterative approximate sparse solver reducing the problem from
O(n^3) to O(n) for n x n matrices.
SymPy's matrix implementation is internally sparse by default and so it can
represent this matrix:
def make(m, n):
L = Matrix(m, n, lambda x, y: -2 if x == y else 1 if x in (y-1,y+1)
else 0)
I = eye(m, n)
kron = KroneckerProduct
# This is not very efficient for sparse case
def kron(A, B):
m, n = A.shape
rows = []
for i in range(m):
row = []
for j in range(n):
row.append(A[i, j]*B)
rows.append(Matrix.hstack(*row))
return Matrix.vstack(*rows)
return kron(I, L) + kron(L, I)
In [12]: %time M = make(100, 100)
CPU times: user 4.47 s, sys: 9.59 ms, total: 4.48 s
Wall time: 4.48 s
In [13]: %time M_squared = M*M
CPU times: user 199 ms, sys: 4.01 ms, total: 204 ms
Wall time: 202 ms
There are only 50000 nonzero entries:
In [14]: M._rep.nnz()
Out[14]: 49600
That makes sense because it is 10000 x 10000 but most rows have only 4 or 5
nonzero elements. This is how the structure of the matrix looks for a
smaller example:
In [19]: make(5, 5).print_nonzero()
[XX X ]
[XXX X ]
[ XXX X ]
[ XXX X ]
[ XX X ]
[X XX X ]
[ X XXX X ]
[ X XXX X ]
[ X XXX X ]
[ X XX X ]
[ X XX X ]
[ X XXX X ]
[ X XXX X ]
[ X XXX X ]
[ X XX X ]
[ X XX X ]
[ X XXX X ]
[ X XXX X ]
[ X XXX X ]
[ X XX X]
[ X XX ]
[ X XXX ]
[ X XXX ]
[ X XXX]
[ X XX]
Whether or not your calculation can be done even with a sparse matrix
implementation depends on what operation you want to perform. The SymPy
implementation can square the matrix using exact integers in 0.2 seconds on
this machine. It could also do the same with arbitrary precision floats but
it is not clear why you want to use bfloats when the numbers are all
integers (you haven't said what you want to do with the matrices).
SciPy's sparse matrix implementation would require using machine precision
int/float types but would otherwise be faster:
In [34]: coords, vals = zip(*make(100, 100).todok().items())
In [35]: i, j = zip(*coords)
In [36]: import scipy.sparse as sp
In [37]: Ms = sp.coo_matrix((vals, (i, j)), (10000, 10000), float)
In [38]: Ms
Out[38]:
<COOrdinate sparse matrix of dtype 'float64'
with 49600 stored elements and shape (10000, 10000)>
In [39]: %time Ms_squared = Ms @ Ms
CPU times: user 18.1 ms, sys: 3.27 ms, total: 21.3 ms
Wall time: 48 ms
Timings vary but it seems to take about 10-20 milliseconds to square the
matrix (at least 10x faster than SymPy). Other operations would likely show
bigger speed differences and SciPy also has a range of different
representations suitable for different operations:
https://docs.scipy.org/doc/scipy/reference/sparse.html
--
Oscar
On Sun, 11 Aug 2024 at 12:53, Barton Willis via Maxima-discuss <
max...@li...> wrote:
> Hi,
>
> Are each of your matrices banded with constant diagonals? If true, you
> could you reduce some questions about your matrices to constant coefficient
> recursion relations.
>
>
> --Barton
> ------------------------------
> *From:* Richard Fateman <fa...@gm...>
> *Sent:* Saturday, August 10, 2024 17:21
> *To:* T Rex <six...@gm...>
> *Cc:* <max...@li...> <
> max...@li...>
> *Subject:* Re: [Maxima-discuss] Error with large matrices
>
> Caution: Non-NU Email
>
> I suggest you try your calculation on a sequence of much smaller matrices
> and see if there is a pattern in your
> computation.
> It seems unlikely that some unique pattern first appears at size 10,000.
>
>
> On Sat, Aug 10, 2024 at 1:53 PM T Rex <six...@gm...> wrote:
>
> Hello everyone,
>
> I'm doing some calculations in which I build a very large square matrix
> full of numbers, I use the "kronecker_product" function to generate it. I
> get an error, to fix it I have to decrease the size of the matrix. The
> size of the matrix should ideally be 10000x10000, when running the RAM
> increases to 1G and the process does not complete. I leave what Maxima
> shows and the code I try to run:
> ////////////////////////////////
> ////////////////////////////////
> Loading ..../.maxima/maxima-init.mac
> Maxima 5.47.0 https://maxima.sourceforge.io
> <https://urldefense.com/v3/__https://maxima.sourceforge.io__;!!PvXuogZ4sRB2p-tU!Bn1euvRjqXJqyJf52rrhzO796mmsp1I4HvTFQ4V2AqROfobjWgkIAkt3Aee6s-ybJYdcrOIH87l7Xg$>
> using Lisp SBCL 2.4.1
> Distributed under the GNU Public License. See the file COPYING.
> Dedicated to the memory of William Schelter.
> The function bug_report() provides bug reporting information.
> (%i1) load(".../ejm.mac");
> Heap exhausted during garbage collection: 0 bytes available, 16 requested.
> Immobile Object Counts
> Gen layout fdefn symbol code Boxed Cons Raw Code SmMix Mixed
> LgRaw LgCode LgMix Waste% Alloc Trig Dirty GCs Mem-age
> 2 0 42 0 135 17 17257 35 0 2 10
> 0 0 0 0.1 566869104 348101034 17321 1 1.3929
> 3 1 67 27 81 14 14598 2 0 1 3
> 0 0 9 0.1 478834032 2000000 14627 0 0.0000
> 4 0 0 0 0 0 0 0 0 0 0
> 0 0 0 0.0 0 2000000 0 0 0.0000
> 5 0 0 0 0 0 0 0 0 0 0
> 0 0 0 0.0 0 2000000 0 0 0.0000
> 6 830 25674 27385 28421 477 184 13 4 34 20
> 25 0 63 2.4 26218128 2000000 11 0 0.0000
> Tot 831 25783 27412 28637 508 32039 50 4 37 33
> 25 0 72 0.2 1071918400 [99.8% of 1073741824 max]
> GC control variables:
> *GC-INHIBIT* = true
> *GC-PENDING* = true
> *STOP-FOR-GC-PENDING* = false
> Collection trigger variables:
> dynamic_space_size = 1073741824
> bytes_allocated = 1071918400
> auto_gc_trigger = 633580067
> bytes_consed_between_gcs = 53687091
> fatal error encountered in SBCL pid 3 tid 3:
> Heap exhausted, game over.
>
> Welcome to LDB, a low-level debugger for the Lisp runtime environment.
> (GC in progress, oldspace=2, newspace=3)
> ldb>
>
> ////////////////////////////////////////////////////////
> ///////////////////////////////////////////////////////CODE
> jLimit:100$
> iLimit:100$
>
> cond1(x,y):= if x=y then ( bfloat(-2) ) else (
> if x=y-1 or x-1=y then (bfloat(1)) else (bfloat(0))
> )$
> L:apply('matrix,makelist(makelist(cond1(i,j),i,1,iLimit),j,1,jLimit))$
>
> cond2(x,y):= if x=y then ( bfloat(1) ) else (bfloat(0))$
> I:apply('matrix,makelist(makelist(cond2(i,j),i,1,iLimit),j,1,jLimit))$
>
> KP_L:kronecker_product(I,L)+kronecker_product(L,I)$
> _______________________________________________
> Maxima-discuss mailing list
> Max...@li...
> https://lists.sourceforge.net/lists/listinfo/maxima-discuss
> <https://urldefense.com/v3/__https://lists.sourceforge.net/lists/listinfo/maxima-discuss__;!!PvXuogZ4sRB2p-tU!Bn1euvRjqXJqyJf52rrhzO796mmsp1I4HvTFQ4V2AqROfobjWgkIAkt3Aee6s-ybJYdcrOKesBOp1g$>
>
> _______________________________________________
> Maxima-discuss mailing list
> Max...@li...
> https://lists.sourceforge.net/lists/listinfo/maxima-discuss
>
|
