nlib: a numerical library written out of a textbook

A library whose documentation is a book

The mdipierro/nlib repository holds two things: the book Annotated Algorithms in Python, subtitled With applications in Physics, Biology, and Finance, and the nlib library that the book builds from the ground up.

That relationship is the whole design. The text says the complete book in PDF is available under a Creative Commons BY-NC-ND License, and it is also available in printed form. The library is described as a pure Python library for numerical computations, and the important constraint is stated plainly: it does not require the module numpy.

So the reader is expected to do two things. Read the derivation in the book, and then read the implementation in nlib, which is the code the chapters are built around. The library is not an incidental companion, it is the worked example.

The licence split follows the same logic. The book is under Creative Commons BY-NC-ND, so it can be shared for non-commercial purposes but not modified or sold. The code is separate: the README states it was created by Massimo Di Pierro in 2016 under a BSDv3 licence.

The repository itself is small, with README.md, a docs/ directory, and src/ as the top level entries, and the library lives in a single file, src/nlib.py. There is no package manifest, no test suite in the top level listing, and no GitHub releases. The last push to the repository was on 2026-04-04.

Linear algebra without numpy

The usage section starts with an import and a matrix, and the matrix is where the design decision shows.

>>> from nlib import *
>>> A = Matrix([[1,2],[4,9]])
>>> print 1/A 
>>> print (A+2)*A
>>> B = Matrix(2,2,lambda i,j: i+j**2)

Two things are worth noticing. The operators are overloaded, so matrix inversion is written as 1/A rather than as a call to a function, and arithmetic on matrices is written the way you would write it on paper. And the matrix constructor is a factory as well as a literal, since B is built as Matrix(2,2,lambda i,j: i+j**2), which is a two by two matrix whose entries are computed by the given function of the row and column index.

That second form is the one that makes the library useful for the kind of problem the book covers. A matrix defined by a generating function is what you want when the matrix represents a discretised differential operator or a covariance structure, and writing that as a nested list comprehension by hand is exactly the tedium the library removes.

The function list shows how much ground the library covers in linear algebra alone. There is Cholesky decomposition, Jacobi_eigenvalues, and the predicates is_almost_symmetric, is_almost_zero, and is_positive_definite, plus condition_number. There are the iterative inverters invert_bicgstab and invert_minimum_residual, and the differentiation machinery: gradient, hessian, and jacobian.

Solvers, optimizers, and how the precision arguments work

The solver and optimizer sections show a consistent naming scheme across the whole numerical core, and the parameters are worth understanding because they are the interface you actually use.

>>> from math import sin
>>> def f(x): return sin(x)-1+x
>>> x0 = solve_newton(f, 0.0, ap=0.01, rp=0.01, ns=100)
>>> print 'f(%s)=%s ~ 0' % (x0, f(x0))

The named parameters are ap for target absolute precision, rp for target relative precision, and ns for maximum number of steps. That trio appears on the optimizers too:

>>> def f(x): return (sin(x)-1+x)**2
>>> x0 = optimize_newton(f, 0.0, ap=0.01, rp=0.01, ns=100)
>>> print 'f(%s)=%s ~ min f' % (x0, f(x0))

The distinction between the two families is the usual one. A solver finds a root of f, so it wants f(x) close to zero, while an optimizer finds a minimum of f, so it wants the derivative close to zero. The library has solve_newton, solve_bisection, solve_secant, solve_fixed_point, and solve_newton_multi for roots, and optimize_newton, optimize_bisection, optimize_secant, optimize_golden_search, plus optimize_newton_multi and a version described as improved, for optimisation.

There is also a differentiation helper, D, in the function list, and the optimiser example uses it to report the derivative at the optimum, which is the standard way to check that an optimiser actually converged rather than merely stopped.

For integration there is integrate, integrate_naive, and integrate_quadrature_naive, and a QuadratureIntegrator class, so the book covers the numerical methods rather than only the closed form ones.

Statistics, and the finance module

The statistics section is compact and shows the library's style of naming.

>>> x = [random.random() for k in range(100)]
>>> print 'μ      =', mean(x)
>>> print 'σ      =', sd(x)
>>> print 'E[X]   =', E(lambda x:x,    x)
>>> print 'E[X²]  =', E(lambda x:x**2, x)
>>> print 'E[X³]  =', E(lambda x:x**3, x)
>>> y = [random.random() for k in range(100)]
>>> print 'corr(x,y) = ', correlation(x,y)
>>> print 'cov(x,y)  = ', covariance(x,y)

The E function is the interesting one. It takes a lambda and a sample and computes the expectation of that lambda over the data, which is a compact way to express the higher moments shown in the example, and it is the kind of primitive that lets you write a distribution without leaving the library.

There is a bootstrap function and a resample function, and confidence_intervals, which together give the resampling side of statistical inference. The variance and sd functions are there, and the function list also shows compute_correlation alongside correlation, which suggests the library separates the raw computation from a convenience wrapper.

The finance module is the one that reaches outside the numeric core. The YStock class takes a ticker and gives you a current snapshot with a price and a market cap, plus a historical() iterator that yields per-day records with a date, an adjusted close, and a log return.

>>> google = YStock('GOOG')
>>> current = google.current()
>>> print current['price']
>>> for day in google.historical():
>>>     print day['date'], day['adjusted_close'], day['log_return']

The log_return field is the finance-specific detail, since returns compound multiplicatively and a log return is what makes them additive. The related classes in the list are Trader, MCEngine for Monte Carlo, MCG for Monte Carlo with a generator, and Markowitz, which is the portfolio construction approach that pairs naturally with the rest.

PersistentDictionary and the concurrency argument

One component in the library is a genuine piece of systems engineering rather than a numerical algorithm, and it comes with an explicit justification.

>>> d = PersistentDictionary(path='test.sqlite')
>>> d['key'] = 'value'
>>> print d['key']
>>> del d['key']

The object works like a drop-in replacement for any normal Python dictionary except that the data is stored in a sqlite database in a file named test.sqlite, so it is still there if you re-start the program. The README compares it to the shelve module, but points out the difference: shelve files cannot safely be accessed by multiple threads or processes unless locked, and locking the entire file is not efficient.

That is a real limitation of shelve and a real reason to use sqlite as the backing store,

The PersistentDictionary is in the function list alongside the graph algorithms, which is an odd place for it, and that is because the book covers a range of topics rather than only numerics. The graph side includes breadth_first_search and depth_first_search, Dijkstra, Kruskal, Prim with a PrimVertex class, DisjointSets, and the string algorithms lcs, needleman_wunsch, encode_huffman and decode_huffman, and make_maze.

There is also a memoize decorator and a memoize_persistent variant, which fits the book's theme: caching is one of the few things that makes the naive versions of the exponential algorithms practical.

Neural networks, plotting, and the numpy exception

The neural network example is a small two-layer network trained on the XOR pattern, and it is the shortest complete example in the documentation.

>>> pat = [[[0,0], [0]], [[0,1], [1]], [[1,0], [1]], [[1,1], [0]]]
>>> n = NeuralNetwork(2, 2, 1)
>>> n.train(pat)
>>> n.test(pat)
[0, 0] -> [0.00...]
[0, 1] -> [0.98...]
[1, 0] -> [0.98...]
[1, 1] -> [-0.00...]

XOR is the standard example because it is the smallest function a single layer cannot represent, so a network that solves it demonstrates the machinery rather than the data. The constructor takes three numbers, here two inputs, two hidden units, and one output.

The plotting section is the one place where the pure Python rule bends, and the README is explicit about the boundary. nlib plotting requires the module matplotlib and numpy for the Canvas object only, and the plots are chainable.

>>> data = [(x0,y0), ...]
>>> Canvas(title='my plot').plot(data, color='red').save('myplot.png')

The available methods are plot, hist, errorbar, and ellipses. The chainable design is what makes it usable for the physics and finance applications the book targets, where an error bar on a measurement is the point rather than an afterthought. The HAVE_MATPLOTLIB constant in the function list is a flag the library exposes so that code can check at runtime whether the optional dependency is present.

That conditional import pattern is the right way to handle an optional plotting backend, and it is worth noting that it is the only numpy or matplotlib requirement in the entire library.

Editorial conclusion

Adopt nlib if you are learning numerical methods and want to read the implementation next to the derivation, or you need a dependency-free numerical core for a constrained environment. Do not adopt it for production numeric work, since a hand-written pure Python implementation is not what you want behind a performance requirement. Verify first that the code runs on your Python version, because the examples in the documentation use the older print statement form.

Frequently asked questions

Official sources

1. Issues
2. mdipierro/nlib on GitHub
3. README