MATLAB: деление матрицы MEX дает результат, отличный от m-файла

Я использовал инструмент кодера MATLAB для создания MEX-версии матричной экспоненциальной функции, которая будет использоваться в другом наборе функций. Проблема в том, что версия MEX дает результаты, отличные от исходного m-файла.

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

Это строка, в которой возникает проблема:

F = (V-U)\(2*U) + I

Где I — единичная матрица размера V и U.

В чем причина несоответствия, когда файл MEX выполняет матричное деление, и как я могу исправить эту проблему? Можно ли переписать эту строку кода без разделения?


matlab matrix equation-solving mex matlab-coder
person G Boggs    schedule 15.10.2014    source источник
comment
Вау.. это действительно странно. Хотя это медленно, попробуйте сделать это: F = inv(V-U)*(2*U) + I. Выполнение <\> по существу берет обратную матрицу левой части и умножает ее на правую. Оператор <\> обычно используется для неквадратных матриц, чтобы найти решение линейной системы уравнений методом наименьших квадратов.   -  person rayryeng    schedule 15.10.2014
comment
@rayryeng Я попробовал это, и, к сожалению, это не работает. Но разве inv(A) не вызывает деление матрицы? inv(A) = A\eye(size(A)) Возможно, есть другой способ сделать это?   -  person G Boggs    schedule 15.10.2014
comment
Не обязательно. inv(A) сама вычисляет обратную матрицу. Существует другой алгоритм вычисления обратной матрицы вместо вызова оператора <\>. Насколько сильно различаются эти два ответа?   -  person rayryeng    schedule 15.10.2014
comment
Что ж, это интересно, потому что числа значительно малы, а некоторые элементы результирующей матрицы точны. Но некоторые из неправильных могут иметь процентное изменение до более чем 100%.   -  person G Boggs    schedule 15.10.2014
comment
Это очень странно. Пробовали ли вы использовать известные значения V и U, а также знать результат F и сравнивать два метода?   -  person rayryeng    schedule 15.10.2014
comment
@rayryeng Я только что сделал, но не уверен, что это сильно поможет. С меньшей матрицей вроде работает (а может просто повезло), но с 10х10 все равно есть НЕКОТОРЫЕ неверные элементы. Хотя разница небольшая, но она все же есть   -  person G Boggs    schedule 15.10.2014
comment
Взгляните на ответ Амро. Я предполагаю, что это может быть ваша версия MATLAB или это может зависеть от операционной системы. Он смог получить оба ответа относительно близко друг к другу (error < 1e-14)   -  person rayryeng    schedule 15.10.2014
comment
@rayryeng, правильно, но проблема не в этом. Проблема в том, что процентное изменение от одной матрицы к другой очень велико, иногда более 100%!   -  person G Boggs    schedule 16.10.2014

Ответы (2)


arrow_upward
4
arrow_downward

У меня нет проблем с созданием кода C из такой операции.

Вот тестовая функция, которую я пробовал:

myfcn.m

function F = myfcn(U,V)
    I = eye(size(U));
    F = (V-U)\(2*U) + I;
end

вот тестовый скрипт, который мы будем использовать для проверки результатов:

test_myfcn.m

U = rand(5);
V = rand(5);
F = myfcn(U,V);

Я начинаю с запуска инструмента генерации кода (ccoder), создаю новый набор проектов для создания MEX-файла, затем добавляю предыдущую функцию myfnc.m в качестве точки входа. Затем я определяю оба типа входных переменных как:

double (:Inf x :Inf)

который задает матрицу MxN неограниченного размера типа double.

Наконец мы можем построить проект. Это производит myfcn_mex.mexw64.

Тестируя как исходную M-функцию, так и сгенерированную MEX-функцию, я получаю практически одинаковый результат (разница в порядке машинного эпсилон):

>> F = myfnc(U,V);
>> FF = myfcn_mex(U,V);
>> norm(F-FF)
ans =
   1.4391e-14
person Amro    schedule 15.10.2014
comment
Проблема не в том, что разница мала, а в том, что процентное отклонение от того, каким должен быть фактический ответ, велико. Это несоответствие вызывает проблему с функцией/программой, в которой я использую этот файл MEX. - person G Boggs; 16.10.2014
comment
@GBoggs: что вы подразумеваете под процентным изменением? это относительная ошибка? Пожалуйста, покажите конкретный пример.. И всегда помните, что числа с плавающей запятой имеют ограничения. Например, в MATLAB (и, возможно, в любом другом языке программирования, использующем IEEE-754): 1.2 - 0.2 - 1 не совпадает с 1.2 - 1 - 0.2. Поскольку сгенерированный код C следует другим этапам выполнения, нет гарантии, что вы получите точную копию результатов исходной функции MATLAB. - person Amro; 16.10.2014
comment
Под этим я подразумеваю следующее: 100*(mex_out - out)./out, это процентное изменение. И для моего применения этой функции это процентное изменение может достигать разницы более 100%, что слишком много. - person G Boggs; 16.10.2014
comment
@GBoggs: это возвращается к проблеме ограничений с плавающей запятой. Ответ, который вы получаете от MATLAB, не более правильный, чем ответ, который вы получаете от сгенерированного кода C. В некотором смысле они оба неверны, потому что, в конце концов, вы не можете представить бесконечную точность в конечном числе битов, и эти ошибки округления накапливаются с течением времени... Более того, эти ошибки растут по-разному в зависимости от порядка в который вы оцениваете шаги алгоритма (код MATLAB не выполняет ту же самую процедуру, что и код C для решения линейной системы) - person Amro; 18.10.2014

arrow_upward
2
arrow_downward

MATLAB внес небольшое изменение в алгоритм EXPM в версии 2014a. Реализации MATLAB Coder отдельные, и соответствующее изменение в алгоритм генерации кода пока не вносилось, поэтому там возможны некоторые неточности. Новая реализация делает (V - U)\(2*U) + I, а старая — (V - U)\(V + U). Они математически эквивалентны, но в целом дают разное поведение округления.

Насколько я знаю, нет систематической разницы в качестве решений линейных систем в MATLAB по сравнению с MATLAB Coder. Основные алгоритмы по существу эквивалентны, а различия в округлении, как ожидается, появятся из различных неясных источников. Остаток может быть меньше для MATLAB или MATLAB Coder в данном случае. Если разница в решениях велика, это говорит о том, что решаемая задача плохо обусловлена. Я могу рассказать об этом подробнее, если хотите, но это описано в каждом учебнике по численному анализу. Можете привести конкретный пример? По крайней мере, вы можете узнать, что там возвращает cond(V - U), когда ваша задача решена в MATLAB?

person Michael Hosea    schedule 20.10.2014
comment
Интересно. Я как-то пропустил момент, когда в OP упоминается экспоненциальная матрица, поэтому я не установил связь с EXPM. Я помню запись в блоге, сделанная недавно Кливом Молером на эту тему. Изменения связаны с этим? - person Amro; 21.10.2014
comment
Я провел быстрый поиск и обнаружил, что Octave реализует тот же аппроксимация масштаба и квадрата Паде, в то время как Python/SciPy ссылается на более новый статья 2009 года Ника Хайама - person Amro; 21.10.2014
comment
Изменение функции MATLAB EXPM было вызвано краевым случаем, возникшим в приложении теории управления. Однако это изменение позволяет EXPM получить правильный ответ на примере, о котором говорил Клив в своем блоге. - person Michael Hosea; 21.10.2014