Оценка параметров с использованием метода наименьших квадратов в Matlab

У меня есть следующие вопросы: Рассмотрим набор уравнений y=ax+b, где я знаю y и x и хочу оценить a и b с помощью метода наименьших квадратов. Предположим, что у нас есть Y=[y1 ; y2] и
A=[x1 1; x2 1], так что Y=A*[a;b]

По методу наименьших квадратов: B=[a;b]=( transpose(A)*A )^-1*transpose(A)*Y

  1. (A'*A) \ A'*Y и A\Y одинаковые?

  2. Какой метод расчета B лучше всего подходит:

    inv( transpose(A)*A ) *transpose(A)*Y

    (transpose(A)*A) \ transpose(A)*Y

    (A'*A) \ A'*Y

    pinv(A)*Y (рассчитать псевдообратную матрицу)

все вышеперечисленное дает немного разные результаты


matlab least-squares
person Alexandros Mel    schedule 25.02.2018    source источник

Ответы (1)


arrow_upward
2
arrow_downward

Прежде чем разрешить ваши сомнения, замечание обязательно. Если вы хотите транспонировать матрицу с помощью сокращенного оператора... вы должны использовать не ', а .'. Первый — это сокращенный оператор для сопряженного транспонирования, а второй — правильный сокращенный оператор для транспонирования. Хотя они обычно дают один и тот же результат, использование первого с матрицами, содержащими комплексные числа, может испортить ваши вычисления.

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

Y = [2; 4];
A = [3 1; 7 1];

Теперь давайте пошагово. Что касается вашего первого ответа, да, две операции эквивалентны с математической точки зрения и дают в основном один и тот же результат:

>> B = A \ Y

B =
    0.5
    0.5

-----------------------------

>> B = inv(A.' * A) * A.' * Y

B =
    0.500000000000001
    0.5

Небольшая разница, которую вы видите, связана с тем, что INV(A) * b менее точна, чем A \ b, и это четко указано даже интерпретатором кода Matlab, если вы наведете курсор на вызов функции inv, за которым следует умножение (которое должно быть помечено оранжевым выделение предупреждения):

Предупреждение

Это частично отвечает и на ваш второй вопрос, но давайте проведем исчерпывающий тест. Я отбросил вычисление, выполненное с использованием inv(A.' * A) * A.' * Y, так как его рекомендуется избегать. Вот так:

tic();
for i = 1:100000
    B = A \ Y;
end
toc();

tic();
for i = 1:100000
    B = pinv(A) * Y;
end
toc();

tic();
for i = 1:100000
    B = (A.' * A) \ A.' * Y;
end
toc();

Это результат бенчмарка:

Elapsed time is 0.187067 seconds.
Elapsed time is 2.987651 seconds.
Elapsed time is 2.173117 seconds.

Учитывая, что три подхода имеют одинаковую точность... первый, безусловно, и, очевидно, самый быстрый.

person Tommaso Belluzzo    schedule 25.02.2018