Решите линейную систему, которая имеет бесконечно много решений с наклонной чертой влево (\) и lsqminnorm. Сравните результаты с помощью 2 норм решений.

Когда бесконечные решения существуют к, каждый из них минимизирует. Команда наклонной черты влево (\) вычисляет одно такое решение, но это решение обычно не минимизирует. Решение, вычисленное lsqminnorm, минимизирует не только norm(A*x-b), но также и norm(x).

Рассмотрите простую линейную систему с одним уравнением и двумя неизвестными. Эта система является недоопределенной, поскольку существует меньше уравнений, чем неизвестные. Решите уравнение с помощью и наклонной черты влево и lsqminnorm.

A = [2 3];
b = 8;
x_a = A\b

x_a = 2×1

         0
    2.6667

x_b = lsqminnorm(A,b)

x_b = 2×1

    1.2308
    1.8462

Эти два метода получают различные решения, потому что наклонная черта влево только стремится минимизировать norm(A*x-b), тогда как lsqminnorm также стремится минимизировать norm(x). Вычислите эти нормы и поместите результаты в таблицу для легкого сравнения.

s1 = {'Backslash'; 'lsqminnorm'};
s2 = {'norm_Ax_minus_b','norm_x'};
T = table([norm(A*x_a-b); norm(A*x_b-b)],[norm(x_a); norm(x_b)],'RowNames',s1,'VariableNames',s2)

T=2×2 table
                  norm_Ax_minus_b    norm_x
                  _______________    ______

    Backslash                0       2.6667
    lsqminnorm      8.8818e-16       2.2188

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

Покажите, как определение допуска к вычислению ранга в lsqminnorm может помочь задать масштаб задачи так, чтобы случайный шум не повреждал решение.

Создайте матрицу низкого ранга ранга 5 и правый вектор стороны b.

rng default % for reproducibility
U = randn(200,5);
V = randn(100,5);
A = U*V';
b = U*randn(5,1) + 1e-4*randn(200,1);

Решите линейную систему с помощью lsqminnorm. Вычислите нормы A*x-b и x, чтобы проверить качество решения.

x = lsqminnorm(A,b);
norm(A*x-b)

ans = 0.0014

norm(x)

ans = 0.1741

Теперь добавьте небольшое количество шума к матричному A и решите линейную систему снова. Шум влияет на вектор решения x линейной системы непропорционально.

Anoise = A + 1e-12*randn(200,100);
xnoise = lsqminnorm(Anoise,b);
norm(Anoise*xnoise - b)

ans = 0.0010

norm(xnoise)

ans = 1.1216e+08

Причина большой разницы в решениях состоит в том, что шум влияет на приближение низкого ранга A. Другими словами, lsqminnorm обрабатывает маленькие значения на диагонали матрицы R в разложении QR A, как являющегося более важным, чем они. Идеально, эти маленькие значения на диагонали R должны быть обработаны как нули.

Постройте график диагональных элементов матрицы R в разложении QR Anoise. Большое количество диагональных элементов находится на порядке 1e-10.

[Q,R,p] = qr(Anoise,0);
semilogy(abs(diag(R)),'o')

Решение этой проблемы состоит в том, чтобы увеличить допуск, используемый lsqminnorm так, чтобы приближение низкого ранга Anoise с ошибкой меньше, чем 1e-8 использовалось в вычислении. Это делает результат намного менее восприимчивым к шуму. Решением с помощью допуска является очень близко к исходному решению x.

xnoise = lsqminnorm(Anoise, b, 1e-8);
norm(Anoise*xnoise - b)

ans = 0.0014

norm(xnoise)

ans = 0.1741

norm(x - xnoise)

ans = 1.0804e-14

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

Создайте 3 3 матрица, которая имеет ранг 2. В этой матрице можно получить третий столбец путем добавления вместе первых двух столбцов.

A = [1 2 3; 4 5 9; 6 7 13]

A = 3×3

     1     2     3
     4     5     9
     6     7    13

Найдите минимальное решение методом наименьших квадратов нормы к проблеме, где равно второму столбцу в. Задайте флаг 'warn' для lsqminnorm, чтобы отобразить предупреждение, если это обнаруживает, что A имеет низкий ранг.

b = A(:,2);
x = lsqminnorm(A,b,'warn')

Warning: Rank deficient, rank = 2, tol =  1.072041e-14.

x = 3×1

   -0.3333
    0.6667
    0.3333