Trefethen & Bau & MATLAB & Julia: Iterative methods
I’m going to wrap up the long-paused MATLAB versus Julia comparison on Trefethen & Bau by chugging through all the lectures on iterative methods in one post.
I’m back to using gists–not thrilled with any of the mechanisms for sharing this stuff.
- Lecture 32 (sparse matrices and simple iterations)
- Lecture 33 (Arnoldi iteration)
- Lecture 34 (Arnoldi eigenvalues)
These are remarkable mainly in that they have such striking similarity in both languages. Aside from square brackets and working around the 1x1/scalar distinction in Julia, little differs besides the syntax of the eigs command.
One frustration, though. I decided to try an interesting alternative to PyPlot in Julia, the Plots package. Actually Plots tries to be a generalization of and alternative route to using PyPlot/matplotlib. I decided to try the PlotlyJS backend instead, however. It makes lovely graphics with very responsive interaction. Since the rendering is in Javascript, I thought it would be perfectly portable, but you can’t see the output in the gist above, even though it should be embedded in the notebook.
I liked using Plots OK; for the most part it’s just different, not better or worse that I could see. I found it awkward to work with subplots. I ended up creating 4 plots individually and then displaying them in a table using another call to plot. I find MATLAB’s setup more convenient. I also could not figure out how to coax a contour plot with a contour at a specified value, which seems like a big lack.
- Lecture 35 (GMRES)
- Lecture 36 (Lanczos and MINRES)
- Lecture 37 (Conjugate gradients)
- Lecture 40 (Preconditioning)
Again the differences are minor. In sparse and iterative methods I found Julia to place a greater emphasis on keyword arguments. For example,
(xCG,~,~,~,resnorm) = cg(A,b,tol=1e-14,maxIter=100);There are default values for tol and maxIter, but if you want to override them you must type the keyword. On the other hand, MATLAB’s arguments are purely positional:
[xCG,~,~,~,resnorm] = pcg(A,b,1e-14,100);If I wanted to specify the maximum number of iterations without changing the default tolerance, then I would need to use an empty matrix in the third position. When one uses a command that does take named parameters as inputs, it’s typically done using 'propname',propval pairs. Except when it isn’t, such as for ODEs and optimization. Confusing! As a user I don’t love typing out the keywords, but Julia at least lets me skip the quote marks. I also know from experience that Julia’s version is a lot easier and clearer to implement on the other side.
So that’s that. I feel that I am at least ready to get off the bunny slopes with Julia. I haven’t found a compelling reason to switch to it, aside from supporting open source software for science (no small thing). Of course I’ve barely scratched the surface. On the flip side, MATLAB has a lot of well-designed and -maintained packages, and its environment still makes a smoother experience for newcomers. If you can afford it, it’s still a great option for interactive numerical computing.
I wonder about the future of Julia. Had Python not gotten a head start, I could see an outpouring of effort to make high-quality Julia packages and Julia being a complete MATLAB reboot. But numpy and scipy do exist, and despite their flaws, they have a huge first-mover advantage. It’s a snap to use Python packages in Julia, so there’s not a dichotomy here. But if the package you want to use a lot exists only in Python, the case for Julia weakens. Overall though, it’s a nice thing that we have several strong, expressive high-level environments for numerical computing. Happy coding!