Anonymous user
Talk:QR decomposition: Difference between revisions
→Questionnable intent
No edit summary |
|||
| (17 intermediate revisions by 6 users not shown) | |||
Line 12:
I believe the code has a problem. Under certain circumstances, division by zero will occur in the line
vdiv(e, vnorm(e, m->m), e, m->m);
The reason seems to be that the loop which creates the Householder matrices loops over the number of rows in the input matrix m (m->m), and attempts to extract a column with the same index as the row index from the matrix z (which has the same dimensions as m). This is fine if m->m <= m->n. For purposes of solving over-determined equation systems, however, we will have m->m > m->n. Consequently, we will try to extract non-existing columns from the matrix z, and what they contain will be indeterminate. Should it happen that all column elements are zero, the division by zero referred to above will take place. This of course renders the whole calculation useless. --[[User:Ernstegon|Ernstegon]] 21:22, 17 November 2012 (UTC)
:: You are right, the code was using the wrong subscript of the matrices. Goes to show that ugly code is more prone to bugs. I hope it's fixed now. --[[User:Ledrug|Ledrug]] 07:33, 19 November 2012 (UTC)
:::Thanks for looking at the problem so quickly. However, I still have doubts. Is the result really a correct QR decomposition? Is not the matrix R supposed to be upper triangular, i.e. all elements below the diagonal should be zero? That is not the case with the results produced with the amended code. When I use it for solving an over-determined system of linear equations, I get results that differ by quite a bit from Matlab - more than could ever be explained by numerical issues. I tried the following simple system in Matlab:<br>
:::A = [1.0 -2.0; 1.0 -1.0; 1.0 0.0; 1.0 2.0; 1.0 3.0]<br>
:::B = [-4.0 -1.0 3.0 6.0 13.0]'<br>
:::X = A\B<br>
:::The result is X = [2.162790698 3.093023256]<br>
:::When I try the amended code in a C test program, using the same input data, I get X = [2.156692 3.108271]. Close, but no cigar.<br>
:::The code as it was before your update produced X = [2.162790 3.093023] - in those cases where it did not crash. --[[User:Ernstegon|Ernstegon]] 16:14, 19 November 2012 (UTC)
::::After a closer look, it seems the two loops which you modified should be:<br>
::::<lang c>for (k = 0; k < m->n; ++k)</lang><br>
::::...not k < m->n - 1<br>
::::With that modification, the code produces correct results in my test program. --[[User:Ernstegon|Ernstegon]] 17:20, 19 November 2012 (UTC)
:::::Eh, ok, maybe it's fixed now. --[[User:Ledrug|Ledrug]] 17:51, 19 November 2012 (UTC)
Am I missing something? I am having trouble compiling the code as is. Is it assuming a certain version/standard of C? --[[User:Cantorg|Cantorg]] 09:51, 18 February 2015 (UTC)
== Common Lisp example ==
Line 29 ⟶ 44:
How about better comments on what each Lisp function needs, does, and returns. In programming terms. Then add it to the task description as the pseudocode to follow? --[[User:Paddy3118|Paddy3118]] 14:08, 23 July 2011 (UTC)
:Why not read the [https://en.wikipedia.org/wiki/QR_decomposition Wikipedia page about QR]? Then if it is not enough, why not open a book on numerical analysis, for example Golub & Van Loan? If a book is asking too much, there are many good online courses about matrix decompositions, [http://www.giyf.com/ Google is your friend]. I'm not claiming that writing a good QR code is easy: it needs some thoughts about geometry, basic matrix operations, and some good ideas about storage (in order to write the Q and R factor in a single n x n matrix plus a vector, as is done in linpack or lapack) -- for example, the Ada solution is not very usable as is, it's just showing how things work with Householder matrix H, in an inefficient way. But programming is not about translating someone else's code, and sometimes it requires some "homework" to understand how to write it yourself. [[User:Arbautjc|Arbautjc]] ([[User talk:Arbautjc|talk]]) 19:55, 6 October 2013 (UTC)
==Draft status==
Line 47 ⟶ 63:
: Are there any reasons to keep this draft? The task appears to be clear enough (mathematical, but that's just its nature) and there are multiple correct implementations. –[[User:Dkf|Donal Fellows]] 13:07, 23 September 2011 (UTC)
:: My only complaint is that it's very difficult to grok, with the way the task description is laid out. But that's something someone can come along and clean up. I agree; there's no substantive reason to keep this as a draft. --[[User:Short Circuit|Michael Mol]] 13:37, 23 September 2011 (UTC)
===Mathematica and Matlab===
The Mathematica and Matlab examples simply show to use the builtin functions for QR-decomposition.
Is this enough? Or should they implement the algorithm given in the description?
I think they should show the algorithm implemented not just the function called!
== Questionnable intent ==
The task is misleading at best. Some answers successfully compute Householder projections,
but, like the task description, they fail to understand that a QR decomposition is never
computed this way, as this would be too much time- and space- consumming.
Only the vector ''u'' of array ''I - s uu''' is ever stored, together with ''s'' or something equivalent to ''s''.
And both ''Q'' and ''R'' are stored in ''A'' in the process, with only a supplementary vector.
LINPACK and LAPACK are of course no exception, though they handle the question in a slightly
different way [http://stackoverflow.com/questions/3031215/mystified-by-qr-q-what-is-an-orthonormal-matrix-in-compact-form].
And when solving a system given an already computed QR in packed form, there are also ways to do
it in a clever way, not computing the ''Q'' and ''R'' matrices effectively.
Most answers only show bad ways to do the job, and they are just showing programming language
in (bad) action.
While I understand that Rosetta Code is not the place for state-of-the-art algorithms, it should be at least
mentionned that they are indeed very poor in this case. And anyway, a reasonnable QR algorithm is not more difficult
to implement, it just requires some work.
[[User:Arbautjc|Arbautjc]] ([[User talk:Arbautjc|talk]]) 11:05, 30 April 2015 (UTC)
| |||