Bug 15211 - PATCH eigen() hangs in R 2.15.2 on both Windows 64-bit and linux-redhat
Summary: PATCH eigen() hangs in R 2.15.2 on both Windows 64-bit and linux-redhat
Alias: None
Product: R
Classification: Unclassified
Component: Low-level (show other bugs)
Version: R 2.15.2
Hardware: x86_64/x64/amd64 (64-bit) All
: P5 major
Assignee: R-core
Depends on:
Reported: 2013-02-17 22:56 UTC by Francois Rousset
Modified: 2013-09-10 10:33 UTC (History)
2 users (show)

See Also:

input matrix (1.03 KB, application/octet-stream)
2013-02-17 22:56 UTC, Francois Rousset
Windows configuration (1.04 KB, text/plain)
2013-02-17 22:57 UTC, Francois Rousset
Linux configuration (905 bytes, text/plain)
2013-02-17 22:57 UTC, Francois Rousset
Patch for the bug (983 bytes, patch)
2013-08-28 20:37 UTC, Aleksey Vorona
Details | Diff

Note You need to log in before you can comment on or make changes to this bug.
Description Francois Rousset 2013-02-17 22:56:17 UTC
Created attachment 1411 [details]
input matrix

eigen() appears to hang indefinitely when calculating eigenvectors for the attached RData file:

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:
Comment 1 Francois Rousset 2013-02-17 22:57:29 UTC
Created attachment 1412 [details]
Windows configuration
Comment 2 Francois Rousset 2013-02-17 22:57:59 UTC
Created attachment 1413 [details]
Linux configuration
Comment 3 Francois Rousset 2013-07-07 19:09:54 UTC
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.

Comment 4 Peter Dalgaard 2013-07-09 10:45:48 UTC
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

          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.
Comment 5 Aleksey Vorona 2013-08-28 20:37:37 UTC
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
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!
Comment 6 Brian Ripley 2013-09-10 10:33:05 UTC
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.]