linear-algebra

Matrix power in R

Trying to compute the power of a matrix in R, I found that package expm implements the operator %^%. So x %^% k computes the k-th power of a matrix. > A<-matrix(c(1,3,0,2,8,4,1,1,1),nrow=3) > A %^% 5 [,1] [,2] [,3] [1,] 6469 18038 2929 [2,] 21837 60902 9889 [3,] 10440 29116 4729 but, to my surprise: > A [,1] [,2] [,3...

Determinant of a complex matrix in R

Is there a way to calculate the determinant of a complex matrix? F4<-matrix(c(1,1,1,1,1,1i,-1,-1i,1,-1,1,-1,1,-1i,-1,1i),nrow=4) det(F4) Error in determinant.matrix(x, logarithm = TRUE, ...) : determinant not currently defined for complex matrices library(Matrix) determinant(Matrix(F4)) Error in Matrix(F4) : complex matrices not...

Large-scale pseudoinverse

I would like to compute the Moore–Penrose pseudoinverse of an enormous matrix. Ideally, I would like to do it on a matrix that has 23 million rows and 1000 columns, but if necessary I can reduce the number of rows to 4 million by only running on one part of my experiment. Obviously, loading the matrix in to memory and running SVD on it...

Matrix Decomposition

I have a square n*n matrix S that has to be decomposed into a product of two matrices - A1 and A2, where A2 is transposed matrix to A1 (A2=A1^T) , so A1 * A2 = S. Are there any algorithms to do such operation effectively? C#/C++ solution would be nice. ...

Solving linear systems of equations

I have a system of 6 equations that I need to solve over and over again in a program (with many different inputs of course). I am currently using the Cramer's rule method of solving the system and it works quite well (it seems that my processor really likes add and multiply operations, it gets solutions in 1 microsecond despite the expli...

Feature selection using Gram-Schmidt orthogonalization in R

Is there any package in R that contains algorithm for feature selection using Gram-Schmidt orthogonalization? ...

Looking for an elegant and efficient C++ matrix library

Greetings, googling for that subject brings, e.g., MTL, exmat, LAPACK and also here. I also seem to remember that Microsoft Research released one, but can't put my hands on it. I look for advice from someone who actually used (or developed...) one of those, hoping to achieve a Matlab experience inside C++ (as much as possible). Thanks in...

Projection of polygon onto plane using GSL in C/C++

The general problem is projecting a polygon onto a plane is widely solved, but I was wondering if anybody could make some suggestions for my particular case. I have a planar polygon P in 3-space and I would like to project it onto the plane through the origin that is orthogonal to the unit vector u. The vertices of P and the coordinates...

Building a 3x3 reflection matrix using GSL

Based on the documents http://www.gnu.org/software/gsl/manual/html_node/Householder-Transformations.html and http://en.wikipedia.org/wiki/Householder_transformation I figured the following code would successfully produce the matrix for reflection in the plane orthogonal to the unit vector normal_vector. gsl_matrix * reflection = gsl...

eigenvectors when A-lx is singular with no solution

How is R able to find eigenvectors for the following matrix? Eigenvalues are 2,2 so eigenvectors require solving solve(matrix(c(0,1,0,0),2,2)) which is singular matrix with no solution. > eigen(matrix(c(2,1,0,2),2,2)) $values [1] 2 2 $vectors [,1] [,2] [1,] 0 4.440892e-16 [2,] 1 -1.000000e+00 > solve(matrix(c(0,1,0,0),...

number of distinct eigenvectors in R

There are 2 eigenvectors corresponding to 1 eigenvalue (with multiplicity 4) for the following example. However, R returns 4 distinct eigenvectors. It looks like pairs of them are approximately the same only differing in machine floating point error (epsilon). Can you please check and confirm? > B [,1] [,2] [,3] [,4] [1,] 2 0 ...

OpenGL Superbible linear algebra - is this correct?

I recently started reading OpenGL Superbible 5th edition and noticed the following: Having just taken linear algebra this seemed odd to me. The column vector is of dimension 4x1 and the matrix is 4x4, how is it possible to multiply them together? If the vector were a row-vector and the output were a row vector I agree that it would ...

normalizing matrices in R

How do I normalize/scale matrices in R by column. For example, when I compute eigenvectors of a matrix, R returns: > eigen(matrix(c(2,-2,-2,5),2,2))$vectors [,1] [,2] [1,] -0.4472136 -0.8944272 [2,] 0.8944272 -0.4472136 // should be normalized to [,1] [,2] [1,] -1 -2 [2,] 2 -1 The function "scale" subtract...

Multivariate Bisection Method

I need an algorithm to perform a 2D bisection method for solving a 2x2 non-linear problem. Example: two equations f(x,y)=0 and g(x,y)=0 which I want to solve somultaneously. I have very familiar with the 1D bisection ( as well as other numerical methods ). Assume I already know the solution lies between the bounds x1 < x < x2 and y1 < y ...

gauss elimination with back substitution code in c

can anyone provide some help with it . i am totally unable to do this one.i tried to do some coding but it didn't work at all. Anyways i have included the code.It is incomplete and wrong anyways. #include<stdio.h> #include<math.h> main() { int rows, columns, iter, i,j,k; printf("enter the number of rows and columns:"); scanf("%d, ...

XNA - About the relation between world space and the screen space

Edit: Just wanted to make the question I have more clear. I pretty much am having trouble seeing how something like Matrix.CreateTransformationZ works in the context of not only matrix multiplication but more importantly what this does to the screen space/world space so I can get a clearer picture. So maybe someone could alter the code o...

Obtaining an invertible square matrix from a non-square matrix of full rank in numpy or matlab

Assume you have an NxM matrix A of full rank, where M>N. If we denote the columns by C_i (with dimensions Nx1), then we can write the matrix as A = [C_1, C_2, ..., C_M] How can you obtain the first linearly independent columns of the original matrix A, so that you can construct a new NxN matrix B that is an invertible matrix with a n...

GSL/BLAS: Multiply a matrix with an inverse matrix

Hi! I'm using the GNU GSL to do some matrix calculations. I'm trying to multiply a matrix B with the inverse of a matrix A. Now I've noticed that the BLAS-part of GSL has a function to do this, but only if A is triangular. Is there a specific reason to this? Also, what would be the fastest way to do this computation? Should I invert A u...

XNA and Matrices

So thanks to many of the replies from this board, I have a much better understanding of some of what's under the hood but I just need this little bit more to get a firm understanding. So I am reading through part of Riemer's tutorials here: http://www.riemers.net/eng/Tutorials/XNA/Csharp/Series2D/Coll_Detection_Matrices.php I'll use th...

[r] LU decomposition with row pivot

The following function does not use row pivoting for LU decomposition. Is there an existing function in R that does LU decomposition with row pivot? > require(Matrix) > expand(lu(matrix(rnorm(16),4,4))) $L 4 x 4 Matrix of class "dtrMatrix" [,1] [,2] [,3] [,4] [1,] 1.00000000 . . ...