Ask Your Question
1

Numerical calculation speed: Sagemath vs. C vs. Matlab

asked 2018-10-10 02:22:46 -0500

DanialBagh gravatar image

updated 2018-10-10 02:25:44 -0500

I mostly use Sagemath for performing symbolic operations on trigonometric polynomials. But I was wondering if it is also a good tool in numerical methods. This includes matrix operations and iterative solution of PDEs with finite difference methods. Usually I take the Sagemath outputs in another program more suited to numerical work. But I was wondering if Sagemath itself is a viable choice in this resepct. Can you tell me about its numerical performance and speed compared to other programs such as C, Matlab, Julia, Mathematica etc.? Thanks

edit retag flag offensive close merge delete

2 answers

Sort by ยป oldest newest most voted
2

answered 2018-10-10 09:08:48 -0500

tmonteil gravatar image

updated 2018-10-10 09:18:18 -0500

I will not answer your question since

  • i am not familiar with non-free Ma* software, so the benchmarks i could do with them will be biased
  • the final choice actually involves a lot of context, like:

    • how fast are you to type unbuggy C code (the implementation time is important too) ?
    • do you plan to apply for jobs that requires you to be fluent in some Ma* ?
    • etc

That said, Sage has numerical capapilities. In particular, it is interfaced with low-level numeric libraries that are as fast as any C code you could write yourself. Somehow, the issue (which is also a strength) is that it also has other capabilities (symbolic, algebraic, ...) so that numerics is not the default, and if you are not careful, then you might work in some other framework and not benefit from numerical software, hence be slower. Here are two advices for numerical computation in Sage.

First advice, the real (or complex) numbers you are working with should explicitely belong to RDF, the Real Double Field, that is, the floating-point numbers that the CPU can handle directly (i mean, not emulated). The default floatin-poing numbers in Sage are emulated with MPFR.

A simpe example:

If you write:

sage: a = 1/10

You define a rational number:

sage: a.parent()
Rational Field

If you write:

sage: a = 0.1

You define a floating-point number:

sage: a.parent()
Real Field with 53 bits of precision

BUT the computations with those numbers (e.g. multiplication) are emulated and do not rely on the CPU floating-point arithmetics. To use such hardware implementation, you need to specify:

sage: a = RDF(0.1)
sage: a.parent()
Real Double Field

Here is a quick benchmark to be explicit, what is going on if you want to elevate a 100x100 matrix to the power 100:

With rational numbers:

sage: d = 100
sage: M = matrix(QQ,d,d)
sage: for i in range(d):
....:     for j in range(d):
....:         M[i,j] = QQ.random_element(distribution='1/n')

sage: %time m=M^100
CPU times: user 15.9 s, sys: 36 ms, total: 15.9 s
Wall time: 15.9 s

This is pretty slow, but actually not ridiculous. Note that the computation uses the blazing fast FLINT library for dealing with rational matrices:

sage: from sage.misc.citation import get_systems
sage: get_systems('M^100')
['FLINT']

However, as it is an exact algebraic computation, the rational numbers that are involved along the computation become very huge, which explains the slowness:

sage: m[42,42]
-143403877757962034908229648979166715246023070908121800149357850500509382879155469212403509993324936983836034231601408025604512188042155054510185506983671055836518268283296097981910299096285130974077147512218469448356315775614176654302917462135918450328466624682433535244279651136503200572437791427649637495142716599330424756533649845742064881459812701278305804645425047594913762093446441972957817257126514503143050407460705950011710084565531827542704863095434609967927914423579742343869924125225487620283312414413274849013393390983736909801421258757308472686338302240842651126642009736715824530821318978456872448162954996631679420804303636154872969963348989833116134054387170749229727283037539004381749906799086536804833060409474136172555366941728618585406059541252034119743871324250678700384872982451171781187582425415036174859973807308791740126700638971051970130984253220733018069676465821436211189459166314816274312856343743847662935629297385039058995867638657628145507773864526108055960240555022481606199202785925058024830432579575255979350088206582498334200765250127346197462183445273244313791481019358399101974261045866228112643969381260915000633692085584764674590962348778981913376088884428003028884525380391147846960296170242506610362212453097147899059246752039470838651568756936447860255647839420934555768524254936307869325709458298713942099571699038352532605173724401046409835891068947995918563622939580020975403667793972445824314924660071421910280960543001126509465348171578165274136712030616591124370278163264756418240035656344485977356981516571170184503490367247329839029918905331047244192646388781483694722070933208276834926440877061903462400310310912714971691777215944382466558195875724235464126574871359191505402431685753780224309872229197968954568292419800681149500774736283850034371556868937848642903100728570736544781094527730083798193064273329641035102107294849482857941069225426202700986705864486343792581539778720880082472939094132176790819879063083321549033366758411669542076025022153052874530561753566578544319391200352505804404263160393897255267328141341573519260317206587567124721064115544033521885878060454490629232621801703261192543213789940147570897598576539236613542839596065958971620169496607381685702579378743276056299101042330609992226050038406371749624774551692227688179116742092296970031968515816849488829812308838973713812548816822625768321587146870870902983652872039699799443046347590284856879915348783897704689183559695047654245372544600457831804001868574295269358336129152534108239647857265940835850541084976890878718168564019735733404413868473483925209291074662970899003410894678407600370270218659911158024027068908727792092608139478540593070795368029820411059560063770818806300201468089597682281464503558441794599205826355228832689993091165952402266711860196750750016774904076511126182712628448362415567227933921041383599916774741258107269767540823396637585176722118472859896814523237395488865515100105081169354315143940698655994219990316698833562289676307780180403165199454157784522590797094478527374783501574320605401121026459046776407888405131418936829305804667907886324195276589175624033831537090441965631102401115345029530651865015732082027908555810312251396237141029959527165808305219771879666805822719914981669182927300885157938700980642128297901433519612881625146340498799586514929278189798256153105789049971719443957612232561467762663752293854016397327300219477211339390143610214020779541658705009870234823049706317349536402840156086539233268937871985432517515463172604709770339182871369732928056025299517733716713747314091327675123334926745606275249259431691897722774610868184641879355032567152931374739472265649186964811897385922336343519610063911873060758568647677747715523847663529768249278160478089091016419362489297291835670381126599721981391398104800068362337473997065285669862259560466340900225586401805671925988517977227517932779620515459839830061885380471763711923114804873916975264575368654385873055716459504494520723831558083395652262993116300740553064319612842833420246247557544382720728988451119281433524116699786263755596727000062127722598041113800958261882375836196955289660885700197622000468013546008085927065493797593246260797703800289284981/3068233910356345990614931617623580052025363483821215700564045905999708872130579435635638399841481262830531624885373508197100019590978701458861914858919475453933335831626337018869221083796738382902937287799220197722602863411596249242017999953591869390022551358931588848413227131713853468406991002164875748927997496445016964500065844454133216282277995902912735256947415434906361273588263278514160727212466551688139678383309809858408034402299238202363975918300188476252228029357782357065266245349761677233243009390356843551637835724941935697834450450391641067171349833358269299691442810496966000919779870369343100391246631460377549774643044365390047051657817585946658234524294649776138187120534351534821840773579829610751310832447486600725060526509550064595122988368195256124864058423717282605793070788109291943826421125164691823812115076026474077869814094507256264903037123942132046194554704463480466194151953222013753575851535909266068351381889442868119162263206899483180358490272017296482144289576513973444294261112059699666422336064124471789391364264348782670434918239459714619940680026205328579847416795834266290123969822991777687122527766743907395330284075498522190911099356582719763293504492190857687113127108868154092846290523404613127281823839319355367841858950096141867125725180537840346174657852219271284649748312949431335428831375228910715855590311463140531395073660981338289678952484137796816476022061933918035644807185484094150372203097845444164474411237557120054395588826596609313324318185619200542474521222601156945937561204053919262758458096622957489337607803141057212019841551394348018575400509414624144762041966797722261413298427247250718642403852941979163565775196630267753594162201644004604128598405799778059296077953215356990298347405719680526218049900875969132764867398447634617212006999958072416302053425176360755489948315972021381329281962092892483968538551199875278932999043228825694147289460821389391542196123732267776705049765257094897877113278229184264926490831631185584733164960076531618470396153473613151234258657456310192170759051863219304456621123355278279853351636831437865030444373485948424006123392441200438124879014148266949522282042492776928390153087373228242848656749639289863488871706330654851583087353991563002244343795515948319346300026716764571683982860817187346244473849237772143213064437969276671050906848978501054587777731162024974665122830861051777737727352881928511637154924639893374845857228747715407851473424933643377720908405529000354330397907209255938805299678934125119669449831465530454086646835933735858140993983339014065247880588078361359731964615456134483883084770134204159814345204244749174449409020049046990520934056010008167850803924410657530397088815782900489352825694993813244668130712579395688762653199029728076253691606927632327043298215051614828668654332776715260887195768963130360672343382240004672189506525524962511197607949891642065867183207363442080766199176149493446122468874516364426863663167075220946179893922974945970201344008638572878764124761956970424348681070095886628516627996606893760461507929409581912516447591129780187151442703213941345129207775426307728902834778231342848802251563985247643798754894859468679778566859578097345891035447137988506394913505228319845793149725530532025618027463700815637740199019814417129437224180845673701989699452489404763133177092415775095597461079206483430742746987539941267231215065034430850591215328608619887106562384176332003631146614379870581183797279760007073870671753318222095559018512362831875252239715243842054200086939822940302835516362615725734811770725841728954988488716481521372103942280119896163860356478135128772559139888825533222699043148399245456222003500672620170732686752410214307258507306619066371241402912392800075101574833897893574114812498440317196603926649564980838400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

Now, with the default floating-point emulated "real" numbers:

sage: M = M.change_ring(RR)
sage: %time m = M^100
CPU times: user 3.63 s, sys: 8 ms, total: 3.64 s
Wall time: 3.64 s

The timing is about 4 times better, but you lose exactness of precision, since the space of representation of numbers stays bounded:

sage: m[42,42]
-4.67382481087635e217

Note that the only external library that is used is MPFR, which is a library for emulated floating-point numbers, not matrices. Hence, the matrix computations use a slow naive Python implementation:

sage: get_systems('M^100')
['MPFR']

Now, let us use the CPU floating-point numbers:

sage: M = M.change_ring(RDF)
sage: %time m = M^100
CPU times: user 4 ms, sys: 0 ns, total: 4 ms
Wall time: 1.08 ms

This is more that 1000 times faster for a very similar result !

sage: m[42,42]
-4.673824810876391e+217

Unfortunately, we do not get any information from the get_systems function:

sage: get_systems('M^100')
[]

This is because there are only Python calls. So, we have to use another tool:

sage: %prun -r m = M^100
         134 function calls in 0.001 seconds

   Ordered by: internal time

   ncalls  tottime  percall  cumtime  percall filename:lineno(function)
        8    0.001    0.000    0.001    0.000 matrix_double_dense.pyx:358(_matrix_times_matrix_)
        8    0.000    0.000    0.000    0.000 matrix_real_double_dense.pyx:90(__cinit__)
        8    0.000    0.000    0.000    0.000 matrix_double_dense.pyx:253(_new)
        1    0.000    0.000    0.001    0.001 power.pyx:105(__pyx_fuse_0generic_power_pos)
        1    0.000    0.000    0.001    0.001 matrix0.pyx:5365(__pow__)
        8    0.000    0.000    0.000    0.000 matrix0.pyx:97(__cinit__)
        8    0.000    0.000    0.000    0.000 matrix1.pyx:2273(matrix_space)
        1    0.000    0.000    0.001    0.001 <string>:1(<module>)
        8    0.000    0.000    0.000    0.000 matrix_double_dense.pyx:116(__create_matrix__)
        8    0.000    0.000    0.001    0.000 element.pyx:3435(__mul__)
        8    0.000    0.000    0.000    0.000 matrix_space.py:1870(nrows)
       16    0.000    0.000    0.000    0.000 matrix0.pyx:4122(is_sparse)
        1    0.000    0.000    0.001    0.001 power.pyx:23(generic_power)
        8    0.000    0.000    0.000    0.000 element.pxd:65(classify_elements)
       16    0.000    0.000    0.000    0.000 matrix_dense.pyx:22(is_sparse_c)
        1    0.000    0.000    0.001    0.001 power.pyx:91(generic_power_long)
        1    0.000    0.000    0.000    0.000 long.pxd:82(integer_check_long)
        1    0.000    0.000    0.000    0.000 long.pxd:195(integer_check_long_py)
        1    0.000    0.000    0.000    0.000 integer.pyx:493(__init__)
        8    0.000    0.000    0.000    0.000 matrix_space.py:1858(ncols)
        8    0.000    0.000    0.000    0.000 element.pxd:108(HAVE_SAME_PARENT)
        1    0.000    0.000    0.000    0.000 integer.pyx:7187(fast_tp_new)
        1    0.000    0.000    0.000    0.000 integer_fake.pxd:48(is_Integer)
        1    0.000    0.000    0.000    0.000 integer.pyx:7265(fast_tp_dealloc)
        1    0.000    0.000    0.000    0.000 integer.pyx:4554(__nonzero__)
        1    0.000    0.000    0.000    0.000 matrix0.pyx:4138(is_square)
        1    0.000    0.000    0.000    0.000 {method 'disable' of '_lsprof.Profiler' objects}
 <pstats.Stats instance at 0x7f4b893ffb00>

Now, we see that the computation relies on _matrix_times_matrix_ metod of matrix_double_dense.pyx, which you can see read online at https://git.sagemath.org/sage.git/tre...

So, we can see that the computation relies on numpy, which is a Python package for numerical matrix operations. This package relies on BLAS, as you can see in the dependencies of numpy: https://git.sagemath.org/sage.git/tre... BLAS is a highly optimized low-level library for dealing with huge floating-point matrices developped by physicists, which is basically what other dedicated numerical software will use.

Second advice, Sage can run any Python package, and there are lots of packages dedicated to numerics. If you want to install a particular Python package to be run within Sage, you can do (from a terminal):

sage --pip install <package>

A pretty exhaustive of Python packages for scientific computation can be found at http://scipy.org/topical-software.html and in particular https://scipy.org/topical-software.ht... for PDEs.

Last, i would be glad if you report any feedback from your experience, so that we could get a better idea and what could be improved in Sage.

edit flag offensive delete link more

Comments

Thank you for the very detailed response. I am trying write a research paper and I want to use free software for doing it, and I think Python based Sagemath Jupyter is gaining traction within the scientific community (did you hear that the economy nobel price winner of this year uses it?). It would benefit me to streamline all parts (both the symbolic computations and the numerical part) in one sheet. Thank you answering the question I was too shy to ask: how to force Sagemath to stop treating numbers as symbolic and just do regular double-precision calculations. I will post the feedback if I did end up writing it in Sagemath instead of taking numerics to a separate C program. Thanks again.

DanialBagh gravatar imageDanialBagh ( 2018-10-10 14:26:53 -0500 )edit
1

answered 2018-10-10 09:12:50 -0500

nbruin gravatar image

For general numerical work, sage includes numpy/scipy. I don't have benchmarks ready for you but these are well-known numerical packages, so with this information you should be able to find comparisons quite easily. If you need to get closer to C, then Cython gives a very good path to write your inner loops with C efficiency with an minimum of interfacing and boiler plate. This is an option not available in matlab. Julia might be able to do that for you. Cython has very good support for numpy's data structures, but some reading will be required to be able to make use of it with full efficiency.

Based on that I would think python (and with that sage) is certainly a viable choice. Whether it's the best one for your application will depend on the expertise available to you (and possibly the strengths of the other packages).

As a matlab alternative there is also octave, which can interface with sage. It will probably depend on which toolkits in matlab you would like to use if that is a viable option.

edit flag offensive delete link more

Comments

Thank you for the answer. I will look up some Numpy benchmarks. I might even write the numerical part in plain python without Sagemath, but I will probably run into problems with regard to lack of linear algebra packages there. But my aim is to distribute the resulting software in a small program. Relying on Sagemath will need the 2-Gig installation.

DanialBagh gravatar imageDanialBagh ( 2018-10-10 14:31:52 -0500 )edit
1

Note that you can show Sagemath program running from a webpage, withour having to ship the whole Sagemath with https://sagecell.sagemath.org/

tmonteil gravatar imagetmonteil ( 2018-10-10 16:48:23 -0500 )edit

Thanks. I think that is meant for very short code (which is a good thing for some situations).

DanialBagh gravatar imageDanialBagh ( 2018-10-11 01:31:21 -0500 )edit

Your Answer

Please start posting anonymously - your entry will be published after you log in or create a new account.

Add Answer

Question Tools

1 follower

Stats

Asked: 2018-10-10 02:22:46 -0500

Seen: 66 times

Last updated: Oct 10