Какова соответствующая функция LAPACK за Matlab [Q, R, E] = qr (A)?

Я в настоящее время пытаюсь дешево вычислить оценку хорошего ранга для матрицы . Поэтому я вычисляю разложение QR-кода с помощью Columnt, используя $A$

[Q,R,E]=qr(A)

в Matlab. Я оцениваю ранг используя $A$

tol = size(A,n)*eps*norm(A,'fro'); 
r = sum(abs(diag(R))>tol)

Это прекрасно работает, и график по всем диагональным элементам R выглядит следующим образом: участок (сортировка (абс (диаг (R)), 1, 'спуск'), 'г *')

Если перенести весь алгоритм на C / Fortran, то я заменяю [Q, R, E] = qr (A), используя DGEQP3 из LAPACK, который также вычисляет столбец, развертывающий декомпозицию QR. Но если я использую ту же оценку для ранга, я в основном получаю что-то не так. Тот же график для созданного из DGEQP3, выглядит так $R$

Матрица ввода одинакова для обоих экспериментов.

Теперь у меня вопрос, на какую функцию LAPACK опирается колонка, поворачивающая QR-разложение от Matlab?

Спасибо за любую помощь, Грису

Изменить: DGEQPF дает тот же неправильный результат.

Edit2:

Входная матрица плотная и построена как $A$ $E+sign(E,F)$
$A$ доступно здесь: http://www-e.uni-magdeburg.de/makoehle/A.mtx.gz (формат MatrixMarket)
Неправильно : http://www-e.uni-magdeburg.de/makoehle/R_wrong.mtx.gz $R$
Я использовал LAPACK 3.4.0 с OpenBlas / GotoBLAS (64 бит)
Matlab 7, 2007b, 2010b Linux 32bit

Edit3: - Используя GDB, я обнаружил, что Matlab 2010b вызывает DGEQP3: # 3 0xaa46ce2f в dgeqp3_ () из /usr/ubuntu10.04/matlabr2010b/bin/glnx86/../../bin/glnx86/../. ./bin/glnx86/mllapack.so Почему я получаю неправильный результат, используя LAPACK (3.4.0 включает исправления, упомянутые в Рабочей записке 176)?

matlab lapack МК ака Грису
источник

Можете ли вы спровоцировать такое же поведение с меньшей матрицей, которой вы могли бы поделиться здесь?

GertVdE

Есть ли у особая структура? Я хочу сказать, что MATLAB использует UMFPACK для разреженной линейной алгебры, но я не уверен, каковы основные библиотеки плотной линейной алгебры.

A

$A$

Арон Ахмадиа

Является ли скудны? Вы можете установить ">> spparms ('spumoni', 1)" и увидеть, что в этом случае matlab использует то, что называется "SuiteSparseQR".

A

$A$

dranxo

Грису - я бы с удовольствием посмотрел на твою матрицу. Однако ссылка www-e.uni-magdeburg.de/makoehle/A.mtx.gz больше не активна (в любом случае, в настоящее время). У вас есть текущая ссылка на матрицу? Спасибо, Les Foster

@LeslieFoster - добро пожаловать в scicomp!

Арон Ахмадиа

Ответы:

Здесь есть две проблемы:

$A$
У вас тот же программный стек, что и во внутренних библиотеках MATLAB?

Плотный или разреженный?

$A$ DGEQP3[Q,R,E] = qr(A) $A$

У вас тот же программный стек, что и во внутренних библиотеках MATLAB?

Вероятно, нет, что может быть одной из причин, почему вы получаете разные результаты.

Я столкнулся с этой проблемой, когда проводил модульное тестирование библиотеки, в которой использовалась QR-факторизация. Я использовал MATLAB для создания прототипа своей работы и получил результаты, отличные от LAPACK или NumPy. Насколько я могу судить, поскольку MathWorks не позволяет легко найти эту информацию, MATLAB использует версию LAPACK не ранее версии 3.1.1 и библиотеку Intel MKL BLAS (для Windows, Intel Mac и Linux) версии 9.1. или выше (см. здесь ). Я не смог ничего найти о версии использования SuiteSparse MATLAB. Покопавшись в Интернете или просмотрев библиотечные файлы для вашей системы, вы сможете найти дополнительную информацию. Вы можете попробовать изменить библиотеки, на которые ссылается MATLAB, чтобы иметь возможность сравнивать с теми же библиотеками в разных пакетах программ; Эрик Чу обеспечивает хорошую рецензиюэто показывает, по крайней мере, как вы можете заменить библиотеку BLAS MATLAB своей собственной (конечно, вы делаете это на свой страх и риск). Он предполагает, что вы можете сделать то же самое с LAPACK. Может даже оказаться возможным заменить версию SuiteSparse, которую использует MATLAB, вашей собственной версией.

$A = QR$ $Q$ $R$

В итоге я использовал NumPy для создания прототипа моих результатов для факторизации QR, потому что он использует системные библиотеки BLAS и LAPACK. NumPy и SciPy не являются заменой MATLAB, так как обеим библиотекам не хватает функциональности MATLAB, но для этой конкретной задачи линейной алгебры Python + NumPy + SciPy + Matplotlib должен работать хорошо.

Джефф Оксберри
источник

Думаю, получить такой же программный стек, как Matlab, невозможно. Использование другой среды для прототипирования также не требуется. Проблема в том, что код работает в Matlab правильно, но не в C.

MK aka Grisu

@Grisu: Я думаю, было бы очень сложно получить такой же программный стек, если не пытаться связать их в своих библиотеках. Что меня смущает, так это то, что вы знаете, что результат в MATLAB верен, а результат в C неверен. Это какая-то тестовая матрица с известными свойствами? Более того, Арон Ахмадиа прав; кроме репликации архитектуры и программного стека, вы не можете ожидать, что получите те же результаты с нестабильным алгоритмом. Мне в основном говорили то же самое на форумах MATLAB два года назад.

Джефф Оксберри

По моему мнению QR не является нестабильным. Я не могу напрямую проверить, какая QR-декомпозиция является правильной, но по рангу и матрице Q я вычисляю проектор, и там я могу легко проверить, хорошие ли результаты получаются хорошими или плохими, а результаты Matlab хороши. Но я пытаюсь связываться с библиотеками Matlab.

МК ака Грису

@ Грису: Есть четкая разница между хорошими результатами и правильными результатами. Недавно я неправильно выполнил расчет, из-за которого мои результаты выглядели великолепно. Тем не менее, расчет был неверным, и правильный расчет сделал мои результаты менее впечатляющими (но, к счастью, показывает, что мои результаты верны). Вы пытаетесь рассчитать ортогональный проектор или наклонный проектор? (Я спрашиваю, потому что значительная часть моей диссертации посвящена наклонным проекторам.)

Джефф Оксберри

@GeoffOxberry: fwiw, на моей версии MATLAB я могу позвонить internal.matlab.language.versionPlugins.blasи internal.matlab.language.versionPlugins.lapackполучить версии BLAS и LAPACK

Amro

См. Страницу Лесли Фостер о программном обеспечении, раскрывающем рейтинг . См. Также эту рабочую записку LAPACK, в которой анализируются сбои в выявлении ранга QRxGEQP3 .

Вы сможете узнать, какие процедуры использует MATLAB, установив точки останова в отладчике и изучив стек. Последнее, что я посмотрел, по общему признанию, несколько лет назад, MATLAB использовал разделяемые библиотеки, и в этом случае имена символов не могут быть удалены, поэтому вы увидите имена функций в стеке вызовов (но не аргументы, потому что они определенно не сохраняют отладочную информацию).

Джед браун
источник

Страница, посвященная программному обеспечению для определения рейтинга, не помогла. Описанная процедура RRQR была первой вещью, которую я использую в своей идее, но она дает еще худшие результаты, чем идея поворота столбцов.

MK aka Grisu

@ Грису - Это должно было помочь тебе. xGEQP3Алгоритм не является полностью безопасным для выявления ранга. Если вы хотите гарантировать, что вы получите правильный результат, вам следует использовать SVD или более безопасный QR-код, такой как xGEQPXили xGEQPY. Нельзя ожидать, что нестабильный алгоритм выдаст одинаковый результат на разных архитектурах или в разных реализациях (возможно, MATLAB использует более старый LAPACK).

Арон Ахмадиа

Я знаю, что GEQP3 не является показателем рангов, но он дает более правильные результаты, чем подпрограммы RRQR. Использование SVD слишком дорого в моем внешнем алгоритме. Я также поговорю с одним из авторов LAWN-176, и он считает, что эта ошибка не устранена. Я также попробовал DGEQPF / DGEQP3 из LAPACK 3.0.0 с теми же результатами.

MK aka Grisu