Talk:LU decomposition: Difference between revisions
m
→Problems in solutions
| (4 intermediate revisions by the same user not shown) | |||
Line 51:
Why is this so?
:First, the LU decomposition is not unique. Notably, the algorithm usually force the main diagonal of L **or** U to be all 1. But there is an additional twist: the "A=LU" decomposition may not exist. This happens if one of the pivots is zero in Gauss algorithm. Then you have to swap rows to continue. And like with Gauss algorithm, it's possible to always apply row pivoting, to get at each step the maximum pivot in absolute value, to enhance the numerical precision. The decomposition is then PA=LU, with P a permutation matrix, thus LU is the same as A with possibly some rows swapped. It's what the VBA code does, as well as other implementations. However, even if the result is on the test cases, there is a bug in the code: the maximum pivot is selected, but not in absolute value. Which means the zero pivot would be selected if all other choices were < 0. I'll change this. [[User:Eoraptor|Eoraptor]] ([[User talk:Eoraptor|talk]]) 13:46, 15 August 2020 (UTC)
::Addition: the VBA code, while not incorrect per se with this modification, is still awfully inefficient (one matrix multiply for each pivot!). All the code needs to be changed. [[User:Eoraptor|Eoraptor]] ([[User talk:Eoraptor|talk]]) 13:51, 15 August 2020 (UTC)
== Problems in solutions ==
After correcting the VBA code, I noticed the Ada, C++ and Lisp codes have the same problem, and probably others as well: when selecting the maximum pivot, it's necessary to get the absolute value, otherwise the zero pivot will be selected if it's the only nonnegative element in the possible choices. This is a **bug**.
There are other problems: the Fortran code did the permutation in a blatantly inefficient way (permute all the lower rows at each step), and the author didn't notice that if you permute the whole rows, you don't have to permute anything afterwards to get the L and U components in A: they are just the lower and upper part of A. It's a matrix trick that's easy to prove.
Overall, I fear several implementations here are badly written. Sadly, it's not the first time I notice this on RC: numerical analysis tasks are not correctly handled. [[User:Eoraptor|Eoraptor]] ([[User talk:Eoraptor|talk]]) 14:21, 15 August 2020 (UTC)
| |||