| Summary: | PATCH eigen() hangs in R 2.15.2 on both Windows 64-bit and linux-redhat | ||
|---|---|---|---|
| Product: | R | Reporter: | Francois Rousset <francois.rousset> |
| Component: | Low-level | Assignee: | R-core <R-core> |
| Status: | CLOSED WONTFIX | ||
| Severity: | major | CC: | pd.mes, voronaam |
| Priority: | P5 | ||
| Version: | R 2.15.2 | ||
| Hardware: | x86_64/x64/amd64 (64-bit) | ||
| OS: | All | ||
| Attachments: |
input matrix
Windows configuration Linux configuration Patch for the bug |
||
Created attachment 1412 [details]
Windows configuration
Created attachment 1413 [details]
Linux configuration
I am not sure why this bug report gets no feedback. The problem has persisted in all versions of R (tested under Windows) up to the current one included, R 3.0.1. F.R. I can reproduce this on OSX Lion with the reference Lapack/BLAS. Apparently not with veclib.
The hanging spot is here, dlapack.f around line 73485 :
BACK = WERR( II )
50 CONTINUE
NEGCNT = DLANEG( N, D, LLD, RIGHT, PIVMIN, R )
IF( NEGCNT.LT.I ) THEN
RIGHT = RIGHT + BACK
BACK = TWO*BACK
GO TO 50
END IF
BACK and RIGHT is Inf and NEGCNT returned as 3 while I is 5.
It is not clear whether this is an issue in R or LAPACK, though. In principle, we are just copying their code, but bugs have crept in before.
A pure FORTRAN version of the issue might be helpful.
Created attachment 1478 [details] Patch for the bug The infinite loop happens in the routine DLARRB(), where estimations are made for a range of values. The loop in question widens the range on each iteration and is unbounded. In the reporter's test case the range grows to INFINITY. The mathematical reason of that is not clear to me. As a result, RIGHT = RIGHT + BACK BACK = TWO*BACK is effectively INFINITY = INFINITY + INFINITY And the loop never ends. I think that the estimation loop can be bounded, but would not risk calculating the proper bound. Instead I have bounded the loop just under IEEE 754-1985 positive infinity (I do not think we can assume ieee_arithmatic is available on all the platforms and it is not the assumption in the rest of LAPACK). I think it is reasonably safe. Using the reporter's test case: > load("/tmp/pb.RData") > bof1 <- eigen(m,only.values=TRUE) > bof2 <- eigen(m) > max(abs(bof1$values - bof2$values)) [1] 5.551115e-16 The computation works. The result is not exactly the same, but the max difference is in the 16th significant digit. I feel that this is a special case and could be detected earlier (not wasting computation cycles going over the estimation loop many many times until it reaches infinity). But my pure mathematics skill is not up to such a task. Hope that helps! Please report LAPACK bug-fixes to LAPACK. Many versions of R use an external LAPACK, and that way they will get the fix too. [You will need to prepare a standalone Fortran version of this example.] |
Created attachment 1411 [details] input matrix eigen() appears to hang indefinitely when calculating eigenvectors for the attached RData file: load("pb.RData") bof <- eigen(m,only.values=TRUE) ## this one works bof <- eigen(m) ## apparent bug bof <- eigen(m[1:90,1:90]) ## this one works bof <- eigen(m[1:91,1:91]) ## apparent bug too I first detected the problem on a linux cluster where a computation that should have taken ~10mn dit not terminate after four days, out of which the example was constructed which reproduces the problem under Windows and linux: Windows configuration: