Fast Algorithms for Finding Pairwise Euclidean Distance (Distance Matrix)

Well, I couldn’t resist playing around. I created a Matlab mex C file called pdistc that implements pairwise Euclidean distance for single and double precision. On my machine using Matlab R2012b and R2015a it’s 20–25% faster than pdist(and the underlying pdistmex helper function) for large inputs (e.g., 60,000-by-300). As has been pointed out, this problem … Read more

Equivalent of `polyfit` for a 2D polynomial in Python

Here is an example showing how you can use numpy.linalg.lstsq for this task: import numpy as np x = np.linspace(0, 1, 20) y = np.linspace(0, 1, 20) X, Y = np.meshgrid(x, y, copy=False) Z = X**2 + Y**2 + np.random.rand(*X.shape)*0.01 X = X.flatten() Y = Y.flatten() A = np.array([X*0+1, X, Y, X**2, X**2*Y, X**2*Y**2, Y**2, … Read more

Solving a linear equation

Cramer’s Rule and Gaussian Elimination are two good, general-purpose algorithms (also see Simultaneous Linear Equations). If you’re looking for code, check out GiNaC, Maxima, and SymbolicC++ (depending on your licensing requirements, of course). EDIT: I know you’re working in C land, but I also have to put in a good word for SymPy (a computer … Read more

Is Google’s Android OpenGL tutorial teaching incorrect linear algebra?

As the guy who wrote that OpenGL tutorial, I can confirm that the example code is incorrect. Specifically, the order of the factors in the shader code should be reversed: ” gl_Position = uMVPMatrix * vPosition;” As to the application of the rotation matrix, the order of the factors should also be reversed so that … Read more

Why does numpy.linalg.solve() offer more precise matrix inversions than numpy.linalg.inv()?

np.linalg.solve(A, b) does not compute the inverse of A. Instead it calls one of the gesv LAPACK routines, which first factorizes A using LU decomposition, then solves for x using forward and backward substitution (see here). np.linalg.inv uses the same method to compute the inverse of A by solving for A-1 in A·A-1 = I … Read more