QR Decomposition with Row Major Matrices

Started by
7 comments, last by Rune Hunter 13 years, 9 months ago
I am attempting to implement QR decomposition with matrices that are row major. The QR decomposition that I am following is on Wikipedia. Doing a direct implantation of the QR decomposition described on Wikipedia yields horribly incorrect results when using row major matrices (i.e. the R matrix is anything but a triangular matrix).

The implementation that works for column major matrices is given below (C# code):

//Orthogonal and normalized matrix of A stored in Q//This method does produce the correct resultsOrthonormalize(ref this, out Q);//Create the R matrix according to Wikipedia//Notice that we are using rows instead of columns because it is row majorR = new Matrix();R.M11 = Vector4.Dot(Q.Row1, Row1);R.M12 = Vector4.Dot(Q.Row1, Row2);R.M13 = Vector4.Dot(Q.Row1, Row3);R.M14 = Vector4.Dot(Q.Row1, Row4);R.M22 = Vector4.Dot(Q.Row2, Row2);R.M23 = Vector4.Dot(Q.Row2, Row3);R.M24 = Vector4.Dot(Q.Row2, Row4);R.M33 = Vector4.Dot(Q.Row3, Row3);R.M34 = Vector4.Dot(Q.Row3, Row4);R.M44 = Vector4.Dot(Q.Row4, Row4);


I have tested this code using Wolfram|Alpha. Because Wolfram|Alpha uses column major matrices (like much of math) I took the input matrix from Wolfram|Alpha, transposed it and sent it through my code. So this code does produce the correct results for column major matrices. The problem is switching the code for row major matrices.

Thanks for any help.
Advertisement
When you say 'row major', are you sure you don't mean row vector? I'm wondering largely because of this:
Quote:Because Wolfram|Alpha uses column major matrices (like much of math)
Generally speaking, linear algebra doesn't know or care about matrix 'majorness'. Many math references do use column vectors, but that's a different issue.

Anyway, getting the terminology cleared up would probably be a good first step towards solving whatever problems you're running into.
Sorry about that. I do mean I am using row vectors while Wolfram|Alpha and the Wikipedia article are using column vectors. As you can probably tell I'm not the greatest at linear algebra yet.
No, Rune Hunter is perfectly fine, row/column major is actually the right term. This is a term used in numeric computing and is due to the fact that different programming languages store matrices in different ways.

This is most commonly encountered when using LAPACK, BLAS (written in FORTRAN) from C, since FORTRAN and C differ in their storage/indexing schemes, which is specified as row/column major format.

I could draw some lousy ASCII diagram here, but instead I can also simply link to the Wikipedia article:
Row/Column major order
Quote:Original post by Hyrcan
No, Rune Hunter is perfectly fine, row/column major is actually the right term. This is a term used in numeric computing and is due to the fact that different programming languages store matrices in different ways.
I know what the terms mean. I was just questioning whether the issue the OP is seeing is really related to matrix layout, or if it's actually due to differing notational conventions. (As far as that question goes, I'm not really arguing one way or another - the OP's post made me wonder, but if you say it's a matrix layout issue, that's good enough for me.)
I did some more messing around after you said that linear algebra often doesn't care or know about row major/column major matrices. When I was performing the orthonomorlization on the matrix to get Q I was using the Gram-Schmidt process. During that process I was using the rows of the matrix as the vectors instead of the columns. After switching that process to use columns I than switched my code I presented above to also use the columns from both Q and A (A being the this pointer). The code now looks like:

//Uses the Gram-Schmidt method using columns instead of rowsOrthonormalize(ref this, out Q);R = new Matrix();R.M11 = Vector4.Dot(Q.Column1, Column1);R.M12 = Vector4.Dot(Q.Column1, Column2);R.M13 = Vector4.Dot(Q.Column1, Column3);R.M14 = Vector4.Dot(Q.Column1, Column4);R.M22 = Vector4.Dot(Q.Column2, Column2);R.M23 = Vector4.Dot(Q.Column2, Column3);R.M24 = Vector4.Dot(Q.Column2, Column4);R.M33 = Vector4.Dot(Q.Column3, Column3);R.M34 = Vector4.Dot(Q.Column3, Column4);R.M44 = Vector4.Dot(Q.Column4, Column4);


Now I correctly get both the Q matrix and the R matrix. The R matrix is now correctly an upper right triangular matrix and when I multiply Q by R I get A. I am not exactly sure why I have to use columns and since my math is being done with row vectors, I'm not sure if there will be any real implications to using columns. When performing orthonormalization (and even orthogonalization) using the Gram-Schmidt process, should I always use the columns of the matrix? Does it matter if the rest of my math uses row vectors?
The standard QR decomposition factors a matrix A as

A = Q1 R1 .

where Q1 is orthonormal and R1 is upper triangular. When you had indexes backwards, are you sure you weren't just factoring it,

A = R2 Q2

with Q2 orthonormal and R2 lower triangular? I.e., (Q2, R2) are the transposes of the QR decomposition of the transpose of A?
Quote:Original post by Hyrcan
No, Rune Hunter is perfectly fine, row/column major is actually the right term.
Hyrcan, I believe you are incorrect. I can't be 100% sure without knowing more about the OP's situation, but I'm pretty sure.

You say that the OP's issue is related to matrix 'majorness', but I propose that it's related to notational convention.

First, note that nowhere in the OP's posted code are 1-d indices used. Instead, we see individual named elements (e.g. M11), and named rows and columns (e.g. Row1, Column1, etc.). Since these member fields identify specific elements, rows, and columns (irrespective of the underlying storage), majorness should not be an issue for the code as shown.

This statement is also telling:
Quote:Original post by Rune Hunter
Because Wolfram|Alpha uses column major matrices (like much of math)
First of all, I can't find anything online that indicates that Mathematica uses column-major matrices; all the references I found suggest that matrices are represented as 'lists of lists', and that they are entered in row-major form by default.

Second of all, it doesn't make much sense to say that 'much of math' uses column-major matrices (in a typical reference on linear algebra, the question of 'row-major vs. column-major' is unlikely to even come up, I would think). However, it is true that much of math makes use of column vectors, another indication that the OP may have been using the wrong term.

Furthermore, the OP himself states that he actually meant 'row/column vector':
Quote:Original post by Rune Hunter
Sorry about that. I do mean I am using row vectors while Wolfram|Alpha and the Wikipedia article are using column vectors.
Now, since there's clearly some confusion about the issue, I suppose the above statement might be incorrect as well, but the evidence doesn't seem to indicate that.

In any case, overall there seems to be more evidence for it being a notation issue (row vs. column vectors) than a matrix layout issue (row- vs. column-major). I don't think there's enough information in the thread to say for sure, but I'd say that my guess is at least as good as yours (and if I had to bet, I'd bet on it being a notational issue).
Jyk, you are absolutely correct that it was a notational issue (if I understand what you mean by that). While implementing the QR decomposition method, I was using Wolfram|Alpha as a reference to see if my method gets the correct answers. It turns out that the main problem is orthonormalizing the matrix. From my understanding there are two ways to orthonormalize a matrix: using column vectors or row vectors. At first I was orthonormalizing the matrix by treating each row of the matrix (M11, M12, etc) as a vector. I would than apply the Gram-Schmidt process to those vectors and reenter the vectors into the matrix as rows again.

As Emergent stated, by treating the rows of the matrix as vectors instead of treating the columns as vectors I basically 'reversed the indices' (for lack of a better name). After a small primitive test I believe I got a solution that was similar to A = R2Q2 where 2 means transpose.

However, now when I orthonormalize the matrix I treat the columns of the matrix as the vectors and, after applying the Gram-Schmidt process to those vectors, I insert the resultant vectors back into the matrix as columns. I than just performed the other operations that Wikipedia gives to get the R matrix and again I treated the columns of the Q and A matrices as vectors. After doing so I was able to correctly obtain the Q and R matrices so that A = QR.

As for the whole Wolfram|Alpha thing, I looked a little more closely and Wolfram|Alpha apparently gives the answer as A = Q^tR where Q^t is the transpose of Q. I guess that really threw me off and I jumped to incorrect conclusions. Also, my lack of knowledge in this specific area didn't help but I'm learning.

Thanks for the help and suggestions everyone!

This topic is closed to new replies.

Advertisement