[patch] [math] changes to QR decomposition
Hi all, I've made some changes to the QR-decomposition code, namely: - rewritten some comments for clarity - extract the Householder transformation routine to a private method - change the decomposition to "economy sized" in which Q is not always square which saves memory in the m >= n case and seems to be how most other libraries handle the decomposition - changed some unit tests Hopefully the extraction of the Householder transformation routine makes subsequent modifications (e.g. adding solver routines) or adaption of the code (e.g. for SVD which was recently discussed?) easier. Joni - To unsubscribe, e-mail: [EMAIL PROTECTED] For additional commands, e-mail: [EMAIL PROTECTED]
Re: [math] Q-R -decomposition
I am reviewing the implementation code and the Householder reflections algorithm to figure out why this is the case. The definitions that I have seen (including the ones that you reference in the javadoc) use the second approach (R is square). Generally, m is assumed to be greater than or equal to n and I think the matrix has to be full rank for the decomp to be unique. I am not an expert in this area, so will defer to others if there are good arguments for using the definition that your implementation uses. The important thing is to document exactly what the code does, both in the interface javadoc and the impl. It would be awkward in this case to have the interface under-specified regarding the dimensions of the factors, so this should probably be settled independently of the algorithm. About the dimensions, it depends on what your requirements are. If you require that the columns of Q form an orthogonal basis of R^m, Q has to be m x m. If that condition is relaxed and m >= n, the last rows of R will be zero. This means that the last columns of Q don't contribute to the product QR. In this case Q can be m x n. The latter form is called "the economy-size decomposition" in MATLAB documentation: http://www.mathworks.com/access/helpdesk/help/techdoc/ref/qr.html The two forms are mentioned also here: http://www.cs.utexas.edu/~plapack/icpp98/node4.html JAMA uses the economy-size decomposition. (After sending my last email I realised that if m <= n, Q can't be m x n already for the fact that at most m vectors can be linearly independent--let alone orthogonal--in R^m.) IIRC, for the factorization to be unique, the matrix has to have full rank and the diagonal elements of R have to be positive, but a factorization exists and can be computed with the Householder algorithm even if the matrix is singular. Another thing I'm not sure of is the naming; I've named the classes without regard to the fact that they work only with real matrices. If there is to be support for complex and arbitrary precision matrices, perhaps a rename would be due? > It seems that many decompositions are used for solving linear systems. > A decomposition object "knows" what the system is like and has access > to the raw factorized data that can be packed in some way, so it would > make some sense to ask it solve systems, too. Plus: users would be > able to switch between different implementations, like with the > Collections API. Beyond what is available in the API (Q and R), what exactly does the QR Decomp "know" that makes solving easier? It has Q stored as a list of Householder vectors. Multiplying a matrix with Q or Q' that is stored in this format is not much more expensive than a normal matrix multiplication, but it's a lot more efficient when compared to first calculating an explicit form of Q and then multiplying. It might be better numerically too, but I cannot be sure.
Re: [math] Q-R -decomposition
> The algorithm used there produces the matrix R and an array of
> Householder vectors. When the getQ() is called, the Householder
> vectors are made into matrices that are multiplied together to yield
> the Q matrix. This seems to be the best way to go about things.
>
That seems fine to me, in terms of the state maintained in the decomp
class and the API as well - i.e., provide the accessors for Q and R,
but maintain state efficiently.
Done; this is what QRDecompositionImpl does now. Singular and
rectangular cases are covered. The tests are more extensive, too.
http://issues.apache.org/jira/browse/MATH-148
There is one thing I'm not sure about: matrix dimensions. Some sources
define the QR-factorization of an m x n matrix so that Q is m x m
(square) and R is m x n, others say that Q is m x n and R is n x n
(square). The current implementation does the former. Of course it's
also possible to define Q as m x r and R as r x n, where r is min{m,
n} or the rank of the original matrix. Do you have any insights on
what should be done?
> I suppose we won't have a base interface for matrix decompositions?
We can talk about that, but if we go with the design above, there is
really no place for it, unless I am missing something. Now is a good
time to discuss these things, though, as once we define the API, we
(and our users) are going to have to live with it for some time.
Other than a "decompose" method (which the design above does not
include), its hard for me to see what a base decomposition interface
would include. The accessors are all different depending on the
decomp. Could be I am missing something, though, so if you have ideas
about how to better structure this, please speak up.
It seems that many decompositions are used for solving linear systems.
A decomposition object "knows" what the system is like and has access
to the raw factorized data that can be packed in some way, so it would
make some sense to ask it solve systems, too. Plus: users would be
able to switch between different implementations, like with the
Collections API.
Then again, solving systems doesn't seem like something inherent to
the idea of a matrix factorization. Could it be a better idea to have
an interface like "LinearSystemSolver" which then can be implemented
by decomposition classes?
Re: [math] Q-R -decomposition
Sorry about the late reply; I was practically all of April on holidays and now I find myself occupied by final exams. +1 - see below. The only real question here is do we need a subpackage for matrix decompositions. Since I think it is unlikely that we will have more than a handful of these, I am OK putting these into the top level, i.e. in .linear. That seems fine to me. > I also had a look at Jama yesterday. There they defer the explicit > generation of the Q part of the decomposition until the user calls > getQ(), which I guess has a computational advantage over calculating > the whole decomp if the user of the API only needs R. This of course > implies that the algorithm has a state and it's most natural to > implement it as a class of its own. Again, I think this should be a separate (immutable) class with state, with the decomp done in the constructor, which should take a RealMatrix (not impl) as argument (using getData to copy if argument is not a RealMatrixImp). I am not sure I understand what you mean about the Q and R accessors in Jama. It looks to me like they are just doing transformations to provide Q and R separately. I think it makes sense to provide those accessors (as we should in the LU case when we externalize that). The algorithm used there produces the matrix R and an array of Householder vectors. When the getQ() is called, the Householder vectors are made into matrices that are multiplied together to yield the Q matrix. This seems to be the best way to go about things. > > From the release plan I read that the QR-decomposition will be needed > for linear regression. Does that mean that it will be used mainly for > least-squares fitting? In that case both Q and R are needed most of > the time, so having the algorithm in a separate class is not strictly > necessary.. The immediate motivation is for solving the normal equations. I don't think we should include the solve() method that Jama has in this class, though. I think it is more natural to have that in the OLS implementation. Tests are a good start. Returning to the overall API design, I think it makes sense to follow the abstract factory pattern used elsewhere in [math] (e.g. the distributions package) to provide for pluggable decomp implementations with defaults provided. So what we would end up with would be an abstract DecompositionFactory class with a concrete DecompositionFactoryImpl subclass providing default implementations. Interfaces for decompositions would be abstracted. User code with look like this: QRDecomposition qr = DecompositionFactory.newInstance().createQRDecomposition(matrix); where QRDecomposition is the interface and DecompositionFactory.newInstance() returns a DecompositionFactoryImpl and createQRDecomposition(matrix) invokes the constructor for QRDecompositionImpl, which is the default implementation. This setup is used in the distributions and analysis packages to provide pluggable implementations. Seems good to me. What do you think is better as a method name, "createQRDecomposition" or "newQRDecomposition"? Both styles seem to be in use. I suppose we won't have a base interface for matrix decompositions? To get started, we can just define QRDecomposition, QRDecompositionImpl. If there are no objections / better ideas, we can then add the factory impls and do the same for LU decomp (and Cholesky, which I think we may also have laying around somewhere). Alright, I'm on it. Joni
Re: [math] Q-R -decomposition
On 3/30/06, Phil Steitz <[EMAIL PROTECTED]> wrote: > Great! The first thing to do is to open a Bugzilla ticket and attach > the code to it, with that apache license in the class file headers > (look at any apache java class for an example). Ideally, you should > also develop and include a test class. Your main method could be the > start of this. Have a look at the RealMatrix test classes for > examples. We can talk further about design here on the list. From a > quick glance at your code, it looks like you have just implemented the > decomp algorithm statically (which is a great start) and we should > talk about how to structure the API, class name, numerical stability, > and package placement. > I have created some tests now and included them in the bugzilla ticket. It would seem most natural to implement QR in math.linear because that's where the rest of the linear algebra related stuff is. Initially I thought the algorithm, as it doesn't require state as such, could be included in the RealMatrix or RealMatrixImpl class, like the LU decomposition. But I'm not sure if that would be pushing too many responsibilities to one class. What is your view on this? I also had a look at Jama yesterday. There they defer the explicit generation of the Q part of the decomposition until the user calls getQ(), which I guess has a computational advantage over calculating the whole decomp if the user of the API only needs R. This of course implies that the algorithm has a state and it's most natural to implement it as a class of its own. >From the release plan I read that the QR-decomposition will be needed for linear regression. Does that mean that it will be used mainly for least-squares fitting? In that case both Q and R are needed most of the time, so having the algorithm in a separate class is not strictly necessary.. - To unsubscribe, e-mail: [EMAIL PROTECTED] For additional commands, e-mail: [EMAIL PROTECTED]
[math] Q-R -decomposition
Hello! I read from the commons math web page that Q-R -decomposition of matrices is on the wish list. Last night and this afternoon I had some free time so I created a clean-room implementation of this decomposition using Householder reflectors. The main source of information was the description of the algorithm found in Wikipedia: http://en.wikipedia.org/wiki/QR_decomposition I would be interested in working towards getting this algorithm integrated in commons-math, but how should I proceed? And does anyone have suggestions on how to improve the implementation? The source code is downloadable from http://www.student.oulu.fi/~jonisalo/QRDecomp.java (should be understandable but may not be up to Apache standards yet) There is a problem with the current implementation though: it only works for non-singular square matrices. What will the usage of the algorithm be like? Should the singular and non-square cases be implemented, too? Kind regards, Joni Salonen - To unsubscribe, e-mail: [EMAIL PROTECTED] For additional commands, e-mail: [EMAIL PROTECTED]
