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
There are quite a few projects that have settled on the Generic Graphics Toolkit for this. The GMTL in there is nice – it’s quite small, very functional, and been used widely enough to be very reliable. OpenSG, VRJuggler, and other projects have all switched to using this instead of their own hand-rolled vertor/matrix math. … Read more
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
This is almost 2x faster for me than your manual version: Reduce(“+”, lapply(names(DT), function(x) DT[[x]] * cf[x])) benchmark(manual = DT[, list(cf[‘A’]*A+cf[‘B’]*B+cf[‘C’]*C+cf[‘D’]*D)], reduce = Reduce(‘+’, lapply(names(DT), function(x) DT[[x]] * cf[x]))) # test replications elapsed relative user.self sys.self user.child sys.child #1 manual 100 1.43 1.744 1.08 0.36 NA NA #2 reduce 100 0.82 1.000 0.58 0.24 NA … Read more
If the matrix A has more rows than columns, then you should use least squares fit. If the matrix A has fewer rows than columns, then you should perform singular value decomposition. Each algorithm does the best it can to give you a solution by using assumptions. Here’s a link that shows how to use … Read more
For a data matrix of size n-by-p, PRINCOMP will return a coefficient matrix of size p-by-p where each column is a principal component expressed using the original dimensions, so in your case you will create an output matrix of size: 1036800*1036800*8 bytes ~ 7.8 TB Consider using PRINCOMP(X,’econ’) to return only the PCs with significant … Read more
Assuming you have two lines of the form Ax + By = C, you can find it pretty easily: float delta = A1 * B2 – A2 * B1; if (delta == 0) throw new ArgumentException(“Lines are parallel”); float x = (B2 * C1 – B1 * C2) / delta; float y = (A1 * … Read more
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
SymPy can calculate arbitrary precision: from sympy import exp, N, S from sympy.matrices import Matrix data = [[S(“-800.21”),S(“-600.00”)],[S(“-600.00”),S(“-1000.48”)]] m = Matrix(data) ex = m.applyfunc(exp).applyfunc(lambda x:N(x, 100)) vecs = ex.eigenvects() print vecs[0][0] # eigen value print vecs[1][0] # eigen value print vecs[0][2] # eigen vect print vecs[1][2] # eigen vect output: -2.650396553004310816338679447269582701529092549943247237903254759946483528035516341807463648841185335e-261 2.650396553004310816338679447269582701529092549943247237903254759946483528035516341807466621962539464e-261 [[-0.9999999999999999999999999999999999999999999999999999999999999999999999999999999999999994391176386872] [ 1]] … Read more
If you are trying to predict one value from the other two, then you should use lstsq with the a argument as your independent variables (plus a column of 1’s to estimate an intercept) and b as your dependent variable. If, on the other hand, you just want to get the best fitting line to … Read more