Tag Info

53

Python (as of 2.6 and 3.0) now searches in the ~/.local directory for local installs, which do not require administrative privileges to install, so you just need to point your installer to that directory. If you have already downloaded the package foo and would like to install it manually, type: cd path/to/foo python setup.py install --user If you are ...

41

I'm going to break up my answer into three parts. Profiling, speeding up the python code via c, and speeding up python via python. It is my view that Python has some of the best tools for looking at what your code's performance is then drilling down to the actual bottle necks. Speeding up code without profiling is about like trying to kill a deer with an ...

32

First, if your undergraduates are like ours and had no prior introduction to computers, expect to spend some time teaching them how to use basic stuff like using a proper editor (i.e., not MS Word), the command line, etc. I think the answer somewhat depends on where you set the focus of your course (or what you are required to teach). For example: How ...

31

This is indeed called catastrophic cancellation. In fact, this particular case is very easy: rewrite the function using the equivalent, numerically stable expression $$\frac{t}{1+\sqrt{1-t^2}}.$$ Since you probably need a reference, this is discussed in most numerical methods textbooks in relation to the formula for solving quadratic equations (that ...

30

Ease of learning Python and Fortran are both relatively easy-to-learn languages. It's probably easier to find good Python learning materials than good Fortran learning materials because Python is used more widely, and Fortran is currently considered a "specialty" language for numerical computing. I believe the transition from Python to Fortran ...

24

I am not super familiar with f2py internals, but I am very familiar with wrapping Fortran. F2py just automates some or all of the things below. You first need to export to C using the iso_c_binding module, as described for example here: http://fortran90.org/src/best-practices.html#interfacing-with-c Disclaimer: I am the main author of the fortran90.org ...

24

Joblib does what you want. The basic usage pattern is: from joblib import Parallel, delayed def myfun(arg): do_stuff return result results = Parallel(n_jobs=-1, verbose=verbosity_level, backend="threading")( map(delayed(myfun), arg_instances)) where arg_instances is list of values for which myfun is computed in parallel. The main ...

22

In 2014, I would've said Python. In 2017, I wholeheartedly believe that the language to teach undergraduates is Julia. Teaching is always about a tradeoff. On one hand, you want to choose something that is simple enough that it is easy to grasp. But secondly, you want to teach something that has staying power, i.e. something that can grow with you. The ...

20

Going from MATLAB to Python does introduce quite a bit of syntax overhead. One way to quantify it is the nice QuantEcon cheatsheet which showcases how there's a lot of extra "stuff" going on when trying to write simple linear algebra commands in Python. The verbose NumPy syntax is really just a symptom of how it was not developed as a technical ...

19

I will try to answer your question considering that you are asking for Python specifically. I will describe my own method of tackling a simulation problem. Strategies for faster simulations are given in this description. First, I prototype new simulations in Python. Of course, I try to make use of NumPy and SciPy as much as I can. Whereas NumPy provides a ...

19

Let me try and break down your requirements: Maintainability Reading/writing text data Strong interfaces/capability for LU factorizations Sparse linear solvers Performance and scalability to large data From this list, I would consider the following languages: C, C++, Fortran, Python, MATLAB, Java Julia is a promising new language, but the community is ...

18

There are two issues that you are likely to be encountering. Ill-conditioning First, the problem is ill-conditioned, but if you only provide a residual, Newton-Krylov is throwing away half your significant digits by finite differencing the residual to get the action of the Jacobian: $$J[x] y \approx \frac{F(x+\epsilon y) - F(x)}{\epsilon}$$ If you ...

17

Here is the Numba solution. On my machine the Numba version is >1000x faster than the python version without the decorator (for a 200x200 matrix, 'k' and 200-length vector 'a'). You can also use the @autojit decorator which adds about 10 microseconds per call so that the same code will work with multiple types. from numba import jit, autojit @jit('f8[:]...

15

A difficulty with any of these types of questions is that the answer is highly community-dependent. To answer some of your questions in haphazard order: MATLAB is used a lot both in academia and in industry. One of the reasons it's used quite a bit in industry is because it is taught in academia. I know for a fact that MATLAB is used at Lincoln Laboratory ...

15

The question has two very different subquestions. I will address the first one only. Matlab's version runs on average 24 times faster than my python equivalent! The second one is subjective. I would say that letting know the user that there is some problem with the integral is a good thing and this SciPy behavior outperforms the Matlab`s one to keep it ...

15

For the first part of my question, I found this very useful comparison for performance of different linear interpolation methods using python libraries: http://nbviewer.ipython.org/github/pierre-haessig/stodynprog/blob/master/stodynprog/linear_interp_benchmark.ipynb Below is list of methods collected so far. Standart interpolation, structured grid: http:/...

15

What you're looking for is Numba, which can auto parallelize a for loop. From their documentation from numba import jit, prange @jit def parallel_sum(A): sum = 0.0 for i in prange(A.shape[0]): sum += A[i] return sum

14

Here is R1, as computed in MATLAB: 1.0e+07 * -7.382605957465515 -9.599867106092937 -2.830412177259742 -0.000000000002830 -0.000000000002830 -1.230434326244253 -1.599977851015490 -0.471735362876624 -0.000000000000472 -0.000000000000472 3.691302978732758 4.799933553046468 1.415206088629871 0.000000000001415 0.000000000001415 -5....

14

First, see Mark L. Stone's answers, which is completely correct. Second, realize that this is the reason why people told you to use relative errors in your numerical analysis class. :) Third, the real question here is why the results do not coincide exactly, since both languages call some BLAS library functions for their computations. There are several very ...

13

Algorithms for rank-1 updates of the SVD (also called incremental SVD) do exist, but I haven't been able to find a LAPACK-like implementation anywhere. The one I've seen mentioned repeatedly is that of Brand (2003). Judging from this website, it seems as though Brand's algorithm is relatively simple to implement using existing LAPACK and BLAS routines as ...

13

I will address only the comparison of C to C++. While it is true that anything written in C can be ported to C++ with a few syntactic touch-ups, the communities have different values. The C library community, more than almost any other, values binary stability. Binary stability is critical for low-level libraries to avoid inflicting constant pain on those ...

13

To the best of my knowledge, Numpy does not support independent streams. Indeed, getting independent streams from the Mersenne Twister (Pythons RNG) is notoriously difficult although it can be done. Consider using the RandomGen package. It is fully compatible with Numpy, and provides you with the PCG64 generator, supporting up to $2^{63}$ independent ...

13

The problem you are encountering is likely not a consequence of your choice of algorithm, but in fact a consequence of the resulting dynamical system after applying time reversal. Per the definition of an attractor, all points in some neighborhood of the attractor will converge to the attractor under the flow of the dynamical system as $t\to\infty$. However, ...

13

There are libraries that you can use in Python that will give you all (or at least nearly all) of the functionality of MATLAB. For example, scipy.integrate.solve_ivp() supports a number of methods for ODE integration that are comparable to what you can get with the various odexxx() functions in MATLAB. So no, you wouldn't have to write your own ODE ...

12

The degeneracy of some eigenvalues looks to me like the hallmark of the breakdown of the Lanczos algorithm. The Lanczos algorithm is one of the more commonly used methods to approximate the eigenvalues and eigenvectors of Hermitian matrices; it's what scipy.eigsh() uses, through a call to the ARPACK library. In exact arithmetic, the Lanczos algorithm ...

11

I wrote a full answer (below the line) before discovering CVXPY, which (like CVX for MATLAB) does all the hard stuff for you and has a very short example almost identical to yours here. You only need to replace the relevant line with p = program(minimize(norm2(A*x-b)),[equals(sum(x),1),geq(x,0)]) My old answer, doing it the harder way with CVXOPT: ...

11

According to the docs, there is no in-place permutation method in numpy, something like ndarray.sort. So your options are (assuming that M is a $N\times N$ matrix and p the permutation vector) implementing your own algorithm in C as an extension module (but in-place algorithms are hard, at least for me!) $N$ memory overhead for i in range(N): M[:,i] = ...

11

CVXOPT only solves (smooth and nonsmooth) convex problems, giving access to several third party convex solvers with guaranteed state of the art worst case complexity. You may pose linear, convex quadratic, linear semidefinite, and many other convex types of constraints. OpenOpt solves general (smooth and nonsmooth) nonlinear programs, including problems ...

11

As Misha and Geoff Oxberry pointed out, Mathematica really has a different focus (just because you can pound in a nail with a screwdriver doesn't mean you should). So I take your question as being "If I know Matlab, why should I learn Python?" [Edit: and so, apparently, did you.] For all intents and purposes, Matlab is the English of scientific computing -- ...

11

Without assuming something special on my_function choosing multiprocessing.Pool().map() is a good guess for parallelizing such simple loops. joblib, dask, mpi computations or numba like proposed in other answers looks not bringing any advantage for such use cases and add useless dependencies (to sum up they are overkill). Using threading as proposed in ...

Only top voted, non community-wiki answers of a minimum length are eligible