Rfast hd.eigen() возвращает NA, но base eigen() не возвращает

У меня проблемы с hd.eigen в Rfast. Он дает очень близкие к eigen результаты для большинства данных, но иногда hd.eign возвращает пустой $vector, NA или другие нежелательные результаты. Например:

> set.seed(123)
> bigm <- matrix(rnorm(2000*2000,mean=0,sd = 3), 2000, 2000)
> 
> e3 = eigen(bigm)
> length(e3$values)
[1] 2000
> length(e3$vectors)
[1] 4000000
> sum(is.na(e3$vectors) == TRUE)
[1] 0
> sum(is.na(e3$vectors) == FALSE)
[1] 4000000
> 
> e4 = hd.eigen(bigm, vectors = TRUE)
> length(e4$values)
[1] 2000
> length(e4$vectors)
[1] 4000000
> sum(is.na(e4$vectors) == TRUE)
[1] 2000
> sum(is.na(e4$vectors) == FALSE)
[1] 3998000

Помимо того факта, что это нарушает мой сценарий, указывают ли эти NA на более глубокую проблему с моими данными? Или hd.eig не может справиться с некоторыми ситуациями, в которых может справиться сток eigen()? Один лучше другого?

Изменить: по предложению Ральфа я проверил свои версии BLAS, и, похоже, R ищет неправильную версию/не в том месте:

~ $ ldd /usr/lib64/R/bin/exec/R
        linux-vdso.so.1 (0x00007ffeec3b9000)
        libR.so => not found
        libRblas.so => not found
        libgomp.so.1 => /usr/lib64/libgomp.so.1 (0x00007feb27ef2000)
        libpthread.so.0 => /usr/lib64/libpthread.so.0 (0x00007feb27ecf000)
        libc.so.6 => /usr/lib64/libc.so.6 (0x00007feb27cdb000)
        /lib64/ld-linux-x86-64.so.2 => /usr/lib64/ld-linux-x86-64.so.2 (0x00007feb27f7b000)

Кроме того, мне неясно, эквивалентен ли openBLAS BLAS, установленному по умолчанию в других дистрибутивах.

> sessionInfo()
R version 3.6.1 (2019-07-05)
Platform: x86_64-generic-linux-gnu (64-bit)
Running under: Clear Linux OS

Matrix products: default
BLAS/LAPACK: /usr/lib64/libopenblas_nehalemp-r0.3.6.so

редактировать 2: я попробовал тот же пример в системе HPC на базе CentOS и не получил никаких NA. Там sessionInfo() раскрывает:

BLAS/LAPACK: /hpc/packages/minerva-centos7/intel/parallel_studio_xe_2019/compilers_and_libraries_2019.0.117/linux/mkl/lib/intel64_lin/libmkl_gf_lp64.so

Редактировать 3: выражение в hd.eign, которое создает NA, равно

vectors <- tcrossprod(y, t(FF) * L^(-0.5))

в частности, L^(-0.5) производит NaN с индексом 2000

> L[2000]
[1] -1.136237e-12

Однако на двух машинах, где NA не возвращаются, L[2000] положителен (хотя и немного отличается, 5.822884e-14 в системе HPC и 3.022511e-12 на моем компьютере с Windows, на котором запущена сборка Microsoft R.

Редактировать 4: разница, по-видимому, возникает в базовой функции eigen(), которая возвращает одно отрицательное значение из матрицы crossprod() xx на проблемной машине, но не два других. Я сохранил объект xx и открыл его между компьютерами, поэтому я знаю, что ввод для eigen() был точно таким же.

Редактировать 5: я углубился на один уровень и обнаружил, что исходное отрицательное значение исходит из этого утверждения в eigen()

    z <- if (!complex.x) 
      .Internal(La_rs(x, only.values))
    else .Internal(La_rs_cmplx(x, only.values))

Редактировать 6: если я сохраню как CSV, а затем снова открою, проблемный компьютер не выдаст отрицательных собственных значений.

> load("/home/james/nfs-cloud/PanosLab/CircRNA/input_to_La_rs.Rdata")
> r <- .Internal(La_rs(as.matrix(x), only.values = FALSE))
> sum(r$values < 0)
[1] 1
> write.csv(x, "test_for_internal.csv", row.names = FALSE)
> x <- read.csv("test_for_internal.csv")
> r <- .Internal(La_rs(as.matrix(x), only.values = FALSE))
> sum(r$values < 0)
[1] 0

Это дает кому-нибудь ключ? Это ошибка?


person Stonecraft    schedule 14.09.2019    source источник
comment
Я не могу воспроизвести ваши результаты, используя R 3.6.1 и Rfast 1.9.5 в Debian Linux.   -  person Ralf Stubner    schedule 16.09.2019
comment
Интересно. Я также использую R 3.6.1 и Rfast 1.9.5, но на Clear Linux. Я попробую на другом компьютере, когда у меня будет возможность.   -  person Stonecraft    schedule 16.09.2019
comment
Вы также можете попробовать использовать один из контейнеров докеров, предоставленных проектом Rocker. Кстати, я использую OpenBLAS для BLAS/LAPCK и gcc версии 9.2.1.   -  person Ralf Stubner    schedule 16.09.2019
comment
Я попробовал это на другой виртуальной машине (Microsoft R для Windows + Rfast). И действительно, hd.eigen() не произвел NA. Интересно, что происходит с моей установкой.   -  person Stonecraft    schedule 16.09.2019
comment
Однако, когда я запускаю тот же код в Ubuntu 18.04 (также R 3.6 и Rfast 1.9.5), я снова получаю NA. Что здесь может происходить?   -  person Stonecraft    schedule 16.09.2019
comment
Какой BLAS/LAPACK вы используете?   -  person Ralf Stubner    schedule 16.09.2019
comment
Давайте продолжим обсуждение в чате.   -  person Stonecraft    schedule 16.09.2019
comment
Согласно sessioninfo(): BLAS/LAPACK: /usr/lib64/libopenblas_nehalemp-r0.3.6.so Однако я думаю, что некоторые ссылки могут быть неработающими, потому что: ldd /usr/lib64/R/bin/exec/R linux-vdso.so.1 (0x00007ffeec3b9000) libR.so => not found libRblas.so => not found   -  person Stonecraft    schedule 17.09.2019
comment
Нередко libRblas.so не устанавливается в стандартное место. Вот почему скрипт, который фактически запускает исполняемый файл R, устанавливает LD_LIBRARY_PATH. В любом случае, у меня нет идей. Одним из способов продвижения вперед будет способ воспроизвести проблему, например. через докер-контейнер.   -  person Ralf Stubner    schedule 17.09.2019
comment
Хорошо спасибо. Поскольку он превратился в другой вопрос, я создал здесь новую тему: stackoverflow.com/questions/57967678/   -  person Stonecraft    schedule 17.09.2019
comment
Хорошо известно, что перебор с флагами оптимизации может сломать численные алгоритмы. Вот почему вам не следует легкомысленно использовать, например, -ffast-math, так как это позволяет компилятору выполнять некоторые преобразования, которые могут нарушить точность. Может быть такая оптимизация была, к сожалению, сделана в Clear Linux. Может быть, вы можете задокументировать матрицу, в которой ClearLinux не может найти собственные векторы, чтобы в конечном итоге это можно было исправить.   -  person Has QUIT--Anony-Mousse    schedule 19.09.2019


Ответы (1)


Функция hd.eigen в Rfast предназначена только для случая, когда n меньше p.

person Mike    schedule 24.09.2019