Multiplying Matrices Without Multiplying

Yes, basically correct. A couple notes/clarifications for other readers:

- The rows a of A are "random," but in the sense of being drawn from some distribution for which we have a training set--not in the sense of, e.g., "assume everything is gaussian" or some other idealized case.

- The training takes a few seconds to maybe a couple minutes, depending on your training set size. But it does have to be done ahead of time--you can't train after receiving the matrices and still get a speedup. Unless maybe the matrices are unrealistically large.

- The error in the approximation is provably not much larger than the error on the training set with high probability (under standard ML assumptions--finite variance, training and test distributions match), but there's no hard limit on it in general. If you know the range of values in your individual vectors, you could obtain a hard limit, although it would be much looser than what you see in practice.

- Correct that many ML applications tolerate large losses in precision. But caveat is that we won't get the right answer to a full 16 bits in each element (this would be ~99.999% accuracy), but instead more like 90-99.9%, depending on the speed-quality tradeoff chosen. So we can't claim it will be sufficiently accurate in all cases

I'm actually not quite sure what you mean by breaking down into two-dimensional operations. We use operations on pairs of vectors, but nothing is assumed to be two-dimensional, and I don't think we suggest imagining anything as a plane? The vectors can be arbitrary and there's no assumption about curvature. Happy to clarify if I can help.

The surface area to volume ratio is just a limit of the shell volume to total volume ratio as the shell thickness goes to zero. So both should asymptotically scale with higher dimensions in the same way.

Please see this thread started by one of the authors, which I would rather see at the top of the page, above my comment:

https://news.ycombinator.com/item?id=28376439

> It is not generally known that to provide seven- decimal accuracy of the sine function, allowing third-order interpolation, only 15 entries are required [*] ...

That got me thinking: is that so? how do you do that?

Luckily, scipy has all sorts of interpolation tools, so one can just experiment a little.

Third-order interpolation -> that means cubic splines. How many points do you need for a cubic spline approximation of the sine function to have 7 decimal place accuracy (in other words the error should be about 0.5e-8? Well, the sine function being anti-symmetric and periodic, you need to only approximate it on the interval [0, pi/2]. It turns out you need 42 points, not including the ends of the interval (where we don't need to tabulate the values of sine; it is simply 0 at 0 and 1 at pi/2).

     import numpy as np
     from scipy.interpolate import CubicSpline
     xs = np.linspace(0,np.pi/2,44,endpoint=True)
     sines = np.sin(xs)
     spline=CubicSpline(xs, sines)
     ts = np.linspace(0,1,10000,endpoint=True)\*np.pi/2
     print("Max Error is ", np.max(np.abs(spline(ts)-np.sin(ts))))
     >> Max Error is  5.02852541828247e-08

But, we can do better. If we tabulate values of sin, then we have for free the values of the cosine at the complementary angles ( cos(x) = sin(pi/2-x). But the derivative of sine is the cosine, so we can use a cubic spline that matches both the values and the derivatives at a number of points. This is called a Cubic Hermite Spline. Turns out we need 23 points.

     from scipy.interpolate import CubicHermiteSpline
     xs = np.linspace(0,np.pi/2,25,endpoint=True)
     sines = np.sin(xs)
     cosines = sines[::-1] #note that we reuse the sines, 
                           #we don't invoke np.cos
     spline=CubicHermiteSpline(xs, sines, dydx=cosines)
     ts = np.linspace(0,1,10000,endpoint=True)*np.pi/2
     print("Max Error is ", np.max(np.abs(spline(ts)-np.sin(ts))))

     >> Max Error is  4.775709550042251e-08

Now that we got to do Hermite interpolation, why stop here? The second derivative of the sine function is the sine with the sign flipped (oh, man, how do you avoid this pun?). So, you don't need to tabulate any new values. In scipy you can use BPoly to get this higher order spline. Now, you only need to tabulate 4 values of the sine function.

     xs = np.linspace(0,np.pi/2,6,endpoint=True)
     sines = np.sin(xs)
     cosines = sines[::-1]
     spline = BPoly.from_derivatives(xs, np.vstack((sines,cosines,-sines)).T)
     ts =np.linspace(0,np.pi/2,10000) 
     print("Max Error is ", np.max(np.abs(spline(ts)- 
                         np.sin(ts))))


     >> Max Error is  2.057940906574629e-08

Well, we can go all the way now. If we only use the values of the sine function and its first 4 derivatives only at the end of the [0,pi/2] interval, we can find a 9th degree polynomial that matches them.

     spline =  BPoly.from_derivatives([0,np.pi/2], 
        [[0,1,0,-1,0], [1,0,-1,0,1]])
     PSpline = PPoly.from_bernstein_basis(spline)
     P = Polynomial(PSpline.c.T[0][::-1])
     print("Max Error is ", np.max(np.abs(P(ts)- 
       np.sin(ts))))

     >>Max Error is  1.700470619869776e-08

And here's that 9th degree polynomial, if you are curious:

x −0.16666666666666785 x^3+8.881784197001252e-15 x^4+0.008331651428065356 x^5+5.167900831715144e-06 x^6−0.000204626778274708 x^7+3.5525642158307225e-06 x^8+1.8946044087406189e-06 x^9