Trapezoidal rule in log-log space

Consider the problem described in this StackOverFlow post. You have a function with certain smoothness properties that are apparent on a log-log plot. This is often accompanied by a large domain of integration. It seems worthwhile to "integrate in logspace", whatever that means. 

This Jupyter notebook probes this question and makes some recommendations.

Randomized SVD

Dimension reduction is an important problem in the era of big data. SVD is a classic method for obtaining low-rank approximations of data.

The standard algorithm (which finds all the singular values) is one of the most expensive matrix decomposition algorithms.

Companies like Facebook or Google deal with huge matrices (big data can be big). Often, they don't care about finding all the singular values - perhaps only the first 10 or 20. They may also not need exquisite precision in the singular values. Good approximations might do just fine.

Fortunately, there are randomized algorithms for finding SVDs which work on a relatively simple logic. One approximates the range of the matrix, by repeatedly multiplying it with random vectors, and works with with those.

The algorithm is fairly simple to implement:

Figure from Erichson et al

In Octave or Matlab, the code can be implemented in about 10 lines.

The resulting truncated-SVD can be a surprisingly good approximation, which can shave multiple orders of magnitude (mileage improves as matrices get bigger) from computation time.

For python, there are decent implementations of randomized SVD in the sklearn package, and the fbpca package from Facebook. This blog post shows some code to call these routines, and provides some benchmarks.

Python: Orthogonal Polynomials and Generalized Gauss Quadrature

A new (to me) python library for easily computing families of orthogonal polynomials.

Getting standard (generalized) Gauss quadrature schemes is extremely simple. For example to get 13 nodes and weights for Gauss-Laguerre integration, correct up to 50 decimal places:

pts,wts = orthopy.schemes.laguerre(13, decimal_places=50)

The numpy Polynomial package provides similar functionality:

pts, wts = numpy.polynomial.laguerre.laggauss(13)

A nice feature (besides arbitrary precision) is that you can derive custom orthogonal polynomials and quadrature rules. All you need to provide is a weight function and domain of the polynomials. From the project webpage:

import orthopy
moments = orthopy.compute_moments(lambda x: x**2, -1, +1, 20)
alpha, beta = orthopy.chebyshev(moments)
points, weights = orthopy.schemes.custom(alpha, beta, decimal_places=30)

This generates a 10-point scheme for integrating functions over the interval [-1, 1], with weight function \(w(x) = x^2\).

Prony Method

Given N equispaced data-points \(F_i = F(t = i \Delta t)\), where \(i = 0, 1, ..., N-1\), the Prony method can be used to fit a sum of m decaying exponenitals: \[F(t) = \sum_{i=1}^{m} a_i e^{b_i t}. \] The 2m unknowns are \(a_i\) and \(b_i\).

In the Prony method, the number of modes in the exponential (m) is pre-specified. There are other methods, which are more general.

Here is a python subprogram which implements the Prony method.
If you have arrays t and F, it can be called as:

a_est, b_est = prony(t, F, m)

Joints from Marginals: Compilation

For convenience, here is a link to the three blogs in this series in one place.

1. A technique for solving the problem in a special case

2. The reason this technique works

3. The corners/edges of this technique, or how it fails for non-Gaussian marginals

Joint from Marginals: non-Gaussian Marginals

In a previous post, I asked the question if the method described here can be used with non-Gaussian distributions.

Let us explore that by considering two independent zero mean, unit variance distributions that are not Gaussian. Let us sample \(x_1\) from a triangular distribution, and \(x_2\) from a uniform distribution.

We consider a triangular distribution with zero mean and unit variance, which is symmetric about zero (spans -sqrt(6) to  +sqrt(6)). Similarly, we consider a symmetric uniform distribution, which spans -sqrt(3) to  +sqrt(3).

Samples from these independent random variables are shown below.

When we use a correlation coefficient of 0.2, and use the previous recipe, we get correlated random variables with zero mean and the same covariance matrix, but ...

... the marginals are not exactly the same!

This is evident when we increase the correlation coefficient to say 0.5.

The sharp edges of the uniform distribution get smoothened out.

Did the method fail?

Not really. If you paid attention, the method is designed to preserve the mean and the covariance matrix (which is does). It doesn't really guarantee the preservation of the marginal distributions. 

Joint from Marginals: Why?

In the previous blog post, we saw a special example in which we were able to sample random variables from a joint 2D-Gaussian distribution from the marginals and the correlation coefficient.

I listed a simple method, which seemed to work like magic. It had two simple steps:

• Cholesky decomposition of the covariance matrix, C(Y)
• Y = LX, where X are independent random variables

The question is, why did the method work?

Note that the covariance matrix of random variables with zero mean and unit standard deviation can be written as, \(C(Y) = E(Y Y')\), where \(E()\) denotes the expected value of a random variable. Thus, we can write the expected value of the Y generated by the method as, \[\begin{align*} E(Y Y') & = E\left(LX (LX)'\right)\\ & = L E(XX') L' \\ & = L I L'\\ & = LL' = C.\end{align*}.\] Here we used the fact that the covariance of X is an identity matrix by design.

Note that this method preserves the covariance matrix (and hence the standard deviation of the marginals).

Does it preserve the mean?

Yes. \(E(Y) = E(LX) = L E(X) = 0.\)

Do the marginals have to be normal for this method to work? Would this work for any distribution (with zero mean, and unit standard deviation)?

We will explore this in a subsequent blog.

Joint Distribution From Marginals

Consider two dependent random variables, \(y_1\) and \(y_2\), with a correlation coefficient \(\rho\).

Suppose you are given the marginal distributions \(\pi(y_1)\) and \(\pi(y_2)\) of the two random variables. Is it possible to construct the joint probability distribution \(\pi(y_1, y_2)\) from the marginals?

In general, the answer is no. There is no unique answer. The marginals are like shadows of a hill from two orthogonal angles. The shadows are not sufficient to specify the full 3D shape (joint distribution) of the hill.

Let us simplify the problem a little, so that we can seek a solution.

Let us assume \(y_1\) and \(y_2\) have zero mean and unit standard deviation. We can always generalize later by shifting (different mean) and scaling (different standard distribution). Let us also stack them into a single random vector \(Y = [y_1, y_2]\).

The covariance matrix of two such random variables is given by, \[C(Y) = \begin{bmatrix} E(y_1 y_1) - \mu_1 \mu_1 & E(y_1 y_2) - \mu_1 \mu_2 \\ E(y_2 y_1) - \mu_2 \mu_1 & E(y_2 y_2) - \mu_2 \mu_2 \end{bmatrix} = \begin{bmatrix} 1 & \rho \\ \rho  & 1 \end{bmatrix},\] where \(\mu\) and \(\sigma\) refer to the mean and standard deviation.

Method

A particular method for sampling from the joint distribution of correlated random variables \(Y\) begins by drawing samples of independent random variables \(X = [x_1, x_2]\) which have the same distribution as the desired marginal distributions.

Note that the covariance matrix in this case is an identity matrix, because the correlation between independent variables is zero  \(C(X) = I\).

Now we recognize that the covariance matrix \(C(Y)\) is symmetric and positive definite. We can use Cholesky decomposition \(C(Y) = LL^T\) to find the lower triangular matrix \(L\).

The recipe then says that we can draw the correlated random variables with the desired marginal distribution by simply setting \(Y = L X\).

Example

Suppose we seek two random variables whose marginals are normal distributions (zero mean, unit standard deviation) with a correlation coefficient 0.2.

The method above asks us to start with independent random variables \(X\) such as those below.

Cholesky decomposition with \(\rho\) = 0.2, gives us,  \[L = \begin{bmatrix} 1 & 0 \\ 0.1  & 0.9797 \end{bmatrix}.\] If we generate \(Y = LX\) using the same data-points used to create the scatterplot above, we get,

It has the same marginal distribution, and a non-zero correlation coefficient as is visible from the figure above.

Logarithms of Negative Numbers

A plot of log(x) looks something like the following:

As x decreases to zero log(x) approaches negative infinity. For negative values of real x, the log function is undefined. For example, consider the following numpy interaction:

>>> import numpy as np
>>> np.log(1)
0.0
>>> np.log(-1)
__main__:1: RuntimeWarning: invalid value encountered in log
nan

If I try to do the same in Octave, I get something different, and interesting.

octave:1> log(1)

ans = 0

octave:2> log(-1)

ans =  0.00000 + 3.14159i

The answer makes sense if we expand the scope of "x" from real to complex. We know Euler's famous identity, \(e^{i \pi} = -1\). Logarithms of negative numbers exist. They just exist in the complex plane, rather than on the real number line.

Octave's answer above just takes the logarithm of both sides of Euler's identity.

We can make python behave similarly by explicitly specifying the complex nature of the argument. So while log(-1) did not work above, the following works just as expected.

>>> np.log(-1+0j)

3.1415926535897931j

For x < 0, if we plot the absolute value of the complex number, then we get a nice symmetric plot for log(x).

Notes:

• In matlab, the command reallog is similar to np.log

A Visual Interpretation of IVPs

Let us  try to explore the structure of an IVP graphically by considering a specific example:\[\dfrac{dy}{dt} = y - y^2, \quad \quad \quad y(0) = 0.1.\] Suppose, we are interested in the solution \(y(t)\) over the domain \(t \in [0, 10]\).

A jupyter notebook, which accompanies this blog, is available here on github.

click to enlarge

Let us consider the 2D domain (y versus t) on which the solution to the IVP is shown as a thick blue line. This is the solution y(t), which satisfies the initial condition, and the differential equation.

We can look at any point (t,y) on this domain, and ask "what is f(y,t) here?"

Note \(dy/dt = f(y,t)\) is the "slope"; it can be positive, negative or zero. The only restriction  is that it cannot point "backwards" in the direction of negative t.

To visualize the slope at each grid point, we can set the horizontal component \(u = 1\), and the vertical component \(v = f(y,t)\), and normalize by dividing each component by \(\sqrt{u^2 + v^2}\).

Therefore, the function f(y,t) completely defines the field of arrows. The streamlines, if you will.

The initial condition \((t_0, y_0)\) tells us where to start in this field; where to "drop the feather" in the river, to be guided and carried away by the streamlines.

The jupyter notebook linked to above, lets you play around with different problems, and different initial conditions. Take a spin.

Curve-Fitting with Python

The curve_fit function from scipy.optimize offers a simple interface to perform unconstrained non-linear least-squares fitting.

It uses Levenberg-Marquardt, and provides a simple interface to the more general least-squares fitting routine leastsq.

Procedure

• Given a bunch of data (\(x_i, f_i\)) and a fitting function \(f(x; a_1, a_2)\), where the \(a_i\) are the parameters to be fit, and \(x\) is the independent variable
• Convert the math function to a python function. The first argument should be the independent variable; the parameters to be fit should follow
• Import the function curve_fit from scipy.optimize
• Call this routine with the function and data. It returns the best-fit and covariance matrix.

Example

Define function:

def f(x, a1, a2):

    f = np.exp(a1*x) * np.sin(a2*x)

    return f

Generate noisy data with a1 = -1 and a2  = 0.5:

xi = np.linspace(0,5.)

fi = f(xi, -1, 0.5) + np.random.normal(0, 0.005, len(xi))

Perform the fit:

from scipy.optimize import curve_fit

popt, pcov = curve_fit(f, xi, fi)

popt = array([-1.00109124,  0.49108962])

Plot:

plt.plot(xi,fi,'o')

plt.plot(xi, f(xi, popt[0], popt[1]))

F2PY Example: Radial Distribution Function

In a recent post, I discussed the use of F2PY for speeding up parts of python code.

The radial distribution function (RDF) or the pair correlation function \(g(r)\) is a statistical mechanical property which quantifies density variation as a function of distance from a reference particle as a function in a system of particles.

It is a fairly fundamental property that can easily be computed from molecular simulations. The Fourier transform of the RDF is experimentally accesible via light scattering.

For the purposes of this demonstration, we are given a bunch of points in a periodic simulation box. Our goal is to find the pair correlation function \(g(r)\), and compare the relative advantage of using f2py to wrap a Fortran 90 subroutine and a vectorized python implementation of the same algorithm.

The details are in this shared iPython notebook on GitHub. You should be able to both (i) see a static html rendition on the linked site, and (ii) download a jupyter notebook for further experimentation.

The bottom line here was that  eventually both python and Fortran gravitated to a \(N^2\) scaling, where \(N\) is the number of particles. Since the Python code exploits vectorization, its performance is not as bad. Indeed, its relative performance in this case seems to improve as the number of particles increase from 500x worse to only 10x worse as the number of particles was increased from 10 to 10000.

F2PY: Interfacing Fortran and Python

Python is notoriously slow with loops. One of the many ways to speed things up is to incorporate subroutines written in C or Fortran into Python.

F2PY is a program that enables us to just that: "provide a connection between Python and Fortran languages." It lets us wrap Fortran code (you'll need a compiler like gfortran) for use in Python.

You can thus salvage legacy Fortran 77 code, or write "new" Fortran 90 code. The latter is much easier to incorporate.

Basic Steps:

1. Write a subroutine file in Fortran 90. Use "intent" modifier to express interface. Avoid allocatable arrays.
2. Use f2py to generate a python interface/signature file ("pyf" file) and a module that can be imported by python. You may have to edit the signature file sometimes.
3. Compile the python module.

Some resources:

1. Robert Johansson has a very nice set of iPython notebooks for Scientific Computing with python. The first part of the lecture "Lecture-6A-Fortran-and-C.ipynb" provides a gentle introduction to F2PY. If you don't have jupyter/iPython notebooks installed on your computer, a PDF of the entire set of notebooks is also available on his github site.

2. A nice example is available here.

3. The official site also provides a cheat-sheet for three different ways to wrap things.

4. Another nice introduction, that goes on to talk about stuff beyond the basic.

I plan on demonstrating this on a practical example shortly.

Block Averaging: Matlab and Python

I'd written about block averaging to estimate error bars from correlated time-series in a couple of blog posts nearly three years ago. Here are the two posts, which explain the motivation and logic behind this technique (post1 and post2).

I wrote programs to carry out this operation in:

• Octave/Matlab
• Python 

The required input is a data-stream or correlated "x" samples. There are optional flags for turning onscreen printing on or off, and to limit the size of the maximum block. The default choice for maximum block size ensures that the datastream is chopped into 4 blocks and more.

The output is a set of three arrays: the block size, the block variance, and the block mean.

The program prints out the mean and variance corresponding to the largest block sizes on screen. This may or may not be the "best choice". A quick look at the plot (by turning the isPlot flag on) will help ascertain this.

Background Material For Students Entering SciComp

Our department, Scientific Computing, is highly interdisciplinary, and we get graduate students from very different backgrounds.

Here is a list of resources, I recommended incoming grad students to look at before they start.

1. Lectures in Differential Equations and Linear Algebra

These video lectures by Gilbert Strang and Cleve Moler present a quick summary of ODEs and major Linear Algebra topics. 

2. An introduction to Matlab

In addition to a compiled language like C++/Fortran/Java, it is useful to know either Matlab or Python

3. Scientific Computing with Python

This contains a set of extremely useful jupyter or iPython notebooks (PDF link if you don't have jupyter installed) which provide a gentle introduction to Python in Scientific Computing. It also includes several advanced topics (parallel programming, incorporating C and Fortran, version control etc.)

The Savitzky-Golay Filter

In 1964, Abraham Savitzky and Marcel Golay published a paper "Smoothing and differentiation of data by simplified least squares procedures" in Analytical Chemistry, which has been heralded as one of the 10 most influential papers in the journal's history.

The Savitzky-Golay filter (SGF) is a digital filter used to smooth noisy data. The basic idea is to chop the dataset into subsets, and then use a low order polynomial to fit successive subsets.

Implementations are available in Octave/Matlab and in recent versions (>0.16) of scipy for python.

Here is a potential use case. I did a Lennard-Jones melt simulation using LAMMPS, and obtained the following pair correlation function g(r) [click to enlarge].

If you look closely, there is a fair amount of noise due to binning.

Let us use SciPy to smooth the noise.

from scipy.signal import savgol_filter
r, gr = np.loadtxt('gr.dat', unpack=True) # read file from disk
gsm = savgol_filter(gr, 15, 4) # smooth it

The second argument (15) has to be an odd number and is the window or subset size, and the third (4) argument is the degree of polynomial to regress. When I plot the smoothed curve:

plt.plot(rFG,grFG,'.')
plt.plot(rFG,gsm1,label='SG')

If you look closely, again:

One can experiment with the window size and the degree of polynomial. In general, a larger window, and a higher degree polynomial make the curve smoother. The figure below shows a window of size 7 and 31 with a degree 4 polynomial.

Intersection of Ellipses

I was recently pondering about the area of intersection of two ellipses. That is, given two ellipses with known shapes (semi-major and minor axes), centers, and relative or absolute orientation, find the area of intersection.

The general problem turns out to be far more complicated than the simpler problem of area of intersection of two circles.

In particular, as Hughes and Chraibi demonstrate (free preprint, paywalled paper), the number of intersections can be 0, 1, 2, 3, or 4. Eberley also has a nice writeup with derivations (pdf).

Luckily, a software implementation (C++) is available here on GitHub.

A Brief History of Runge-Kutta Methods

Butcher's "A History of Runge-Kutta Methods" contains a beautiful account of the early days of this now ubiquitous method for solving ordinary differential equations.

Carl Runge, a German mathematician and spectroscopist, toyed with the idea of improving upon Euler's simple time stepping approximation for solving ODEs. In a 1895 paper, he wrote three down three different "Runge-Kutta" schemes to solve a particular differential equation that arose in the study of atomic spectra.

He used so much mathematics in his research that physicists thought he was a mathematician, and he did so much physics that mathematicians thought he was a physicist.

Two other German mathematicians, Karl Heun and Martin Kutta made additional contributions. In 1900, Heun wrote down RK methods up to order 4. In 1901, Kutta extended it to order 5, and systematized most of the prior work. His contribution included the classical 4th order RK method, which is the ODE analog of Simpson's method for numerical integration.

Numerical Recipes has this to say about this century old algorithm:

For many scientific users, fourth-order Runge-Kutta is not just the first word on ODE integrators, but the last word as well. In fact, you can get pretty far on this old workhorse, especially if you combine it with an adaptive stepsize algorithm.

By 1925, Nystrom had extended and cleaned up most of the RK methods developed previously.

Over the years, higher order methods were developed, stability properties were studied, and implicit RK methods capable of solving stiff problems were developed.

Warren Buffett 90-10 Portfolio and the Kelly Criterion

When Warren Buffett was asked how he would ask his wife to manage their money after he was gone, his advice was simple.

My advice to the trustee could not be more simple: Put 10% of the cash in short-term government bonds and 90% in a very low-cost S&P 500 index fund. (I suggest Vanguard’s.) I believe the trust’s long-term results from this policy will be superior to those attained by most investors…

Can we try and understand this choice through the lens of the Kelly criterion, which we have explored in some depth in previous posts?

Setup

Suppose that returns \(r\) from the stock market obey a Gaussian distribution, with the mean \(\mu = 1.1\) and scale factor \(\sigma = 0.3\). Note that this \[r = 1 + \text{annual rate of return}.\] This choice essentially assumes average annual returns of 10% from the market, and anticipates typical swings (\(2\sigma\)) to lie between -50% and +70%.

We will play with some of these parameters later, but right now, they seem reasonable enough; so lets not fixate on them.

I will also assume that the worst possible return is \(r = 0\), which implies a return of -100%, which means everything is completely wiped out. This has never happened in the past, and it seems like a good bet that this will "never" happen in the future. Similarly, I will assume that the best possible returns are 300%. It seems like a reasonable cap.

So my annual returns are expected to be random variables drawn from a truncated normal distribution. Let's see what this would look like.

Kelly Criterion

In the Kelly criterion, the crucial expression to develop is the mean of the logarithm of the expected return.

In this case, the probability of obtaining a return of \(r\) is \[\mathcal{N}_{T}(r|\mu, \sigma, r_{min}, r_{max}) dr,\] where the subscript "T" reinforces the "truncated" aspect of the normal distribution.

If I wager a fraction \(f\) of my investment pot, then I should expect to end up with, \(1 - f + f r\), where I assume that the uninvested fraction is kept safe in cash, and earns 0%.

Thus I can write the logarithm of the pay off as:
\[\log S = \langle \log R \rangle = \int_{r_{min}}^{r_{max}}  \mathcal{N}_{T}(r) \log(1 + f(r-1)) dr\]  Let's visualize log S for this choice of \(\mathcal{N}_{T}(r)\).

Python codes used to generate these (and subsequent results) are at the end of this post.

The optimal allocation to the S&P portfolio using the Kelly criterion (by maximizing log S) is 82.8%, which is reasonably close to Buffett's suggested portfolio.

Scenarios

Let us consider three scenarios: (i) conservative (\(\mu = 1.05\), \(\sigma = 0.5\)), (ii) moderate (\(\mu = 1.08\), \(\sigma = 0.4\)), (iii) aggressive (\(\mu = 1.12\), \(\sigma = 0.3\)).

The optimal allocations for the conservative, moderate, and aggressive assumptions are approximately, 55%, 83%, and 100%, respectively.

If we assume a \(\sigma = 0.3\), then a 90% stock allocation corresponds to a \(\mu \approx 1.09\), or an expected 9% rate of return from the S&P.

Historically, the standard deviation of S&P returns has been closer to 0.2. I like to use a higher number because the histogram of the actual returns is not perfectly normal, and the tails have been fatter than one would expect from a \(\sigma = 0.2\).

Python Code

Thoughts on Truncation and Round-off Errors

In scientific computing or numerical analysis, truncation and roundoff errors pop up all over the place. In some operations the former dominates, in others it is the latter, and in a few others, they both conspire in surprising ways.

Truncation error is "is the error made by truncating an infinite sum and approximating it by a finite sum". Or as I sometimes like to say, "somewhere a Taylor series is truncated, and it is not happy!" 

For example, in numerical quadrature, different Newton-Cotes formulas may be derived by locally "fitting" the integrand function with polynomial approximations (Taylor series), and computing the area under the polynomial, since it is much easier to do.

In the solution of ordinary differential equations using Runge-Kutta or linear multistep methods, one tries to perform a term-by-term comparison with a Taylor series expansion, and use polynomial interpolation of the recent past, respectively, to develop different formulas.

In these operations, from a practical standpoint, truncation error dominates!

Roundoff error is the difference between a mathematical number and its finite precision computer representation. It is purely an artifact of the fixed word size used in modern computers. It often dominates the conversation, when subtraction of nearly equal numbers is carried out numerically.

For instance, in solving a system of linear equations, the idea is to take linear combinations of equations to subtract and zero out certain coefficients. In linear systems, there is clearly no Taylor series, and roundoff error is the root cause of ill-conditioning. In fact, if we had infinitely precise number representation on computers, it might be argued that we don't have to worry about conditioning at all!

Since linear systems are at the root of nearly all problems in scientific computing, round-off error is always lurking around the corner.

There are some spectacular cases like numerical differentiation where both roundoff error and truncation error play an "equally" important role. In deriving finite difference formulae, we truncate appropriate Taylor expansions. Furthermore, the derivative itself involves ratios of differences of nearly identical terms.