Teaching Scientific Computing with Julia

Community Teaching & Outreach
1

This post is to continue a discussion posted by @Bruno_Amorim in another thread.

I occasionally teach a class about numerical methods. On one occasion (a small advanced course) I managed to use Julia with my students, but in the other class (a massive 2nd-year course) I would need to first convince my department (I think it will be possible soon enough).

Now, I’ve started an online resource that aims to provide a reference for this matter, which is my Numerical Computing in Julia tutorial booklet.

I made that tutorial in hurry to cover the immediate needs of that class, but I’m wondering of how to continue developing it.

I have various topics on my to-do list: optimization (with AD), finite differences, finite elements, Interval Arithmetic, and so on.

However, besides covering more content, I was wondering if there is some discussion we could have about the packages that we perhaps should consider “standard” in this area, and that we should refer to beginners/students.

Now, I’m not a big fan of massive packages where you define a problem and then you call solve. I’d rather rely on a slightly more modular approach, which also makes the math more transparent, and where access to the actual code being run is easy. Especially for this goal of teaching a class.

For example, for quadrature, I’d rather say in my class that the best option is probably FastGaussQuadrature which provides nodes and weights in O(n) time (this gives something to discuss in the class, as well), and if not applicable (due to unknown singularities), Quadgk is probably the go-to package (but maybe this just reflects my previous Matlab experience and there is a better package?).

For iterative solvers, we also had some discussion on this thread about the differences between KrylovKit, Krylov, and IterativeSolvers. I think KrylovKit has the best generic API, but still couldn’t find the time to analyze it more thoroughly (like, how easy is to include preconditioners). I’m inclined to refer to KrylovKit as the cleanest solution.

As for eigensolvers, I’ve successfully used the eigs function from the LinearAlgebra package. I’m not aware of its limitations. Why should we refer to any other package?

For Plotting, I think Plots.jl is the standard as it has so many backends.

For piecewise polynomial interpolation (splines or piecewise-linear) I think the Interpolations.jl is perfectly fine. For high-order polynomial interpolation, I’d rather share some code snippets and then refer to ApproxFun.

Anyway, that was off the top of my head, any discussion is welcome!

10 Likes
2

Just use LinearSolve.jl. It covers all of this, and more.

Those are all one dimensional. Integrals handles higher dimensional. What we need to do is add a FGQ method in Integrals.jl.

DataInterpolations.jl has the most complete set here. And it has regression splines as well.

It’s not iterative, and doesn’t expose some of the bits for resolving.

It’s fine enough. Though these days I’m starting to lean more and more towards Makie, I wouldn’t teach it yet.

13 Likes
3

Thanks @ChrisRackauckas, it seems I’ve got quite a number of packages to check out

Oh, and this thread also appeared at the same time I wrote this:

https://discourse.julialang.org/t/julia-for-education-package-selection/

so I guess we can continue the discussion there, as well.

3 Likes
4

I must say that I like your choice of topics and problems so far.

However, there is one aspect that is far too rarely the content of lectures and could fit very well into your course. The numerical algorithms you have covered so far can actually all be implemented with a few lines of code. For small problems, this is also completely valid and I would do exactly the same. However, in scientific programming, sooner or later you always reach the point where problems are too complex to be programmed down straightforward and with ease. At this point, whether you have learned it or not (I didn’t in my studies), you have to deal with software engineering and approach problems in a more structured way. For example, extending the functionality of existing types, creating your own types and integrating them into the type hierarchy, using and designing interfaces, writing tests and documentation, and so on. I think some of these aspects could be outlined very well in your lecture. Such skills can also be transferred well when the students are confronted with other programming languages at some point.

We teach a scientific programming course at the Hamburg University of Technology in Germany that has a much stronger focus on software engineering. We try to build a bridge between numerics and software engineering and get a lot of positive feedback from our students.

8 Likes
5

That’s a fantastic approach.

I had to learn software engineering the hard way: making mistakes, throwing away and starting projects over, and trying to get valuable insights from software engineering books that are aimed at commercial software developers, so a lot of time is spent discussing client specifications and so on, which very likely don’t apply for scientific software.

I wrote a very general discussion about software engineering for scientific code in my post, on why it’s both desirable and difficult to learn, and what should be our main concerns. I absolutely will continue writing about this. I think Julia makes it easier to grasp good software engineering practices, thanks to its built-in testing framework and automated documentation, so it’s a good idea to put emphasis on these topics when using Julia.

6 Likes
6

I teach a similar course. I try not to dwell too much on the software engineering aspect (that is covered by other courses) but I (try to…) emphasize the “split the problem in parts, make each part a black box that doesn’t know about the rest, test them separately on simpler problems and then make them talk to each other” part which is important in scientific computing and that doesn’t really come naturally to students (who tend to code everything together, and then are lost when it doesn’t work).

The overarching project is a nonlinear PDE (I use Turing’s model of pattern formation, to change a bit from the mostly mechanics-type problems they have encountered before) solved by an implicit method, where you define an integration method (backwards Euler), a nonlinear solver (Newton), a linear solver, the functions defining the problem, test them separately and then make them talk to each other. Then for extra credit you can swap the simple methods with more efficient black box solvers (diffeq, nlsolve, iterativesolvers…). That’s the theory at least; in practice it proves a bit too ambitious (in particular the fact that the solution has two components that depend on time and space makes the students lose a lot of time converting back and forth between different array formats), and I’m happy if I can at least make them sweat a bit on the debugging and learn the hard way that you have to test parts separately. This is a topic that is a bit of a challenge to teach because any advice basically falls on deaf ears when people haven’t experienced first-hand the challenges, and that takes time.

12 Likes
7

I would say that some sort of adaptive quadrature (e.g. QuadGK or HCubature, perhaps via via Integrals.jl) should probably be the go-to “default” method — usually, you know your error tolerance, not the corresponding number of nodes, and it’s also a safer default to use a method that’s robust to badly behaved integrands in a small portion of the domain.

5 Likes
8

Your blog post reflects exactly my personal experience. I only came into contact with programming during my doctoral thesis. At the beginning Python, then later Julia. I also have to say that there was a barrier in Python that separated me, as a scientist without C++ or C knowledge, from the developers of real packages. I don’t have that feeling in Julia, because everyone speaks the same language.

2 Likes
9

Notice that the function eigs has been moved to the package Arpack.jl. The problems I have faced with eigs occur when using shift-and-inverse. Sometimes you get errors with the factorization step that is used to solve the linear problems. And there is no way to change the factorization used.

ArnoldiMethod.jl and KrylovKit.jl do not implement shift-and-inverse, but it is not easy to adapt them to do so. See for example Spectral transformations section in the docs of ArnoldiMethod.jl. However, it would be nice to have something ready out-of-the-box.

Between ArnoldiMethod.jl and KrylovKit.jl, from experience with a MSc student (computing some bandstructures from a tight-binding Hamiltonian), we found that KrylovKit.jl was more accurate (we would get “noisy” bands when using shift-and-inverse with ArnoldiMethod.jl, for the eigenvalues further away from the target). The downside of KrylovKit.jl is that it does not support in-place operations.

As I said in another post I think that for such a high level operation, such as plotting, being capable of switching backends is not such a great advantage. If you really want to fine-tune the look of a plot (and you always end up wanting to do so), you end up running into the issue of not all backends having the same features, and having to use backend dependent options. At that point, it is just better to use directly the backend.

On the other hand, for other kinds of problem, such as eigenproblems or linear problems, I think that being capable of easily changing backends/solvers is actually a useful thing. You might be writing code to solve a general problem where the eigen/linear-problem is just one of the steps. However, there might not be an ideal solver. It might depend on the size of the problem you are dealing, or you might want to easily spawn the solver for testing. So being capable of passing a solver as a function argument really makes things easy. In this respect, I think that Floops.jl does this really well.

I am not fully convinced that transforming

x = A\b

into

prob = LinearProblem(A, b)
sol = solve(prob)
x = sol.u

Is a great idea for beginners, but I guess one can write it as:

x = solve(LinearProblem(A, b)).u

which is a one-liner, but still a bit verbose.
For me a notation like:

x = linearsolve(A, b, method)

seems more natural (before Julia I mostly used Mathematica). But I guess that when the method needs to use a cache, maybe the best option is really to use something like problem-init-solve.

10

I would actually say the reverse : it is easy enough to get it working, as that tutorial shows. It’s not a one liner, but it’s something you can copy paste into a four line function, by using an orthogonal combination of different packages (krylovkit, a factorization, and linear maps) which you can then swap out (if you need a different factorization for example). Otoh, arpack’s direct support of shift inverse results in a super messy api (i never can figure out what the meaning of the flags is in shift inverse mode). It’s true that it is slightly less convenient to have to write that function, but on the other hand it makes it explicit that a factorization is needed, with all the performance tradeoffs.

11

Not a fan of this, that means you are encoding the type of the problem in the name of the function. It’s similar to calling functions stringprint and floatprint instead of just print and then dispatching on the input type.

Being able to pass a problem object around, and then dispatching on its type seems more ‘Julian’.

(Accessing internal object fields like sol.u feels uncomfortable though.)

4 Likes
12

I do not share that opinion. Yes it is true that it is simple to implement shift-and-inverse by oneself. But then the advantage of having a common interface for eigenproblems for which we can change the method is lost. The use case I have in mind is the following:

I want to compute a banstructure** (a series of eigenvalues that depend on a parameter k, for several values of said parameter). So I would define a function

bandstructue(h, kpts, eigensolver)

where h returns a matrix for each value k stored in kpts. The eigensolver argument determines the method I will use to perform the diagonalization. Very often, I only want inner eigenvalues. So having a shift-and-inverse eigensolver is convenient. I think one can have an API that allows to change the factorization that is used to solve the shift-and-inverses.

The FEAST eigensolver algorithm also requires the solution of linear problems (either with factorizations or iterative solutions in IFEAST). The Jacobi-Davidson algorithm also requires the solution of linear problems. So I think we need to come up with a good way of passing linear solver to eigensolvers.

** for me this comes from a tight-binding Hamiltonian.

13

I guess this is a matter of aesthetic preference. The name solve to me seems too generic to be meaningful. While things such as linearsolve, eigen immediately say what you are doing. Of course we can write solve(LinearProblem(A, b)) which then clearly states what you are doing in one line. So I guess this is fine

Perhaps it could be possible to write:

x, _ = solve(LinearProblem(A, b))

such that it unpacks the solution type? Similar to writing

vals, vecs = eigen(A)

which unpacks the eigenvalues and eigenvectors fields from the Eigen object that is returned by eigen?

14

I don’t see the issue? Just define a function my_eigensolver(A) that does the shift invert transformation and calls krylovkit or some such, and then pass it to bandstructure? (I happen to have the exact same problem as you but with dft hamiltonians, so no shift invert but lazy operators and preconditioners, and that’s what I do : DFTK.jl/diag.jl at c4f03a0f75ecdb5f923f9777df1f40c4648e84c3 · JuliaMolSim/DFTK.jl · GitHub)

15

I agree. I find it more important, that the package has an interface that has a Julian feel to it and this is something Plots offers.

1 Like
16

Just take a LinearSolve type :sweat_smile: .

While you do give up a little bit of simplicity, what you gain in uniformity greatly outweighs it. Iterative solvers aren’t weird, they are just another choice of linear solver. Nonlinear problems aren’t so different from linear problems. The idea of tolerances is the same, the way you control verbosity, etc. all is the same. If you know one you know them all. \ is nice until it’s not.

2 Likes
17

The purpose of Plots.jl’s staying power is the recipe system. Any discussion of Plots and a replacement has to discuss the alternative to handling arbitrary types without recipes.

1 Like
18

That may be true in a scientific context. However, in a lecture where plotting is more of a tool rather than the focus this discussion might not be necessary. At least I never used this feature in an educational context. Nothing I want students to deal with, if it is their first lecture, where they have contact with Julia.

19

Wait what? The whole purpose is for simplification and lectures. Recipes are the solution that makes it so you don’t have to do teach the entire type setup. Why would you want to teach students about the entire type trees of packages in order to make a plot? Why would plot(convert(Vector{Float64},map(real,x)) be better pedagogically than plot(x) allowing the recipes to do the conversion?

7 Likes
20

Those are good points indeed.

I guess in lectures we should keep a focus on simplicity, as well.

Plots, HCubature, LinearSolve, and DataInterpolations, have simple APIs, and seem very useful to be used in lectures.

For eigenvalues, we could say that LinearAlgebra handles full-matrix methods. I still have to re-check all the alternatives for matrix-free methods.