Документация

Вычислительные факторы

Одной из самых важных проблем в техническом вычислении является решение систем одновременных линейных уравнений.

В матричном обозначении общая проблема принимает следующую форму: Учитывая две матрицы A и b, действительно там существует уникальный матричный x, так, чтобы A x = b или x A = b?

Это поучительно, чтобы рассмотреть пример 1 на 1. Например, делает уравнение

имеет уникальное решение?

Ответ, конечно, является да. Уравнение имеет уникальное решение x = 3. Решение легко получено делением:

Решение обычно не получается путем вычисления инверсии 7, который является ^7–1 = 0.142857..., и затем умножение ^7–1 21. Это было бы, больше работают и, если ^7–1 представлен конечному числу цифр, менее точных. Подобные факторы применяются к наборам линейных уравнений с больше чем одним неизвестным; MATLAB^® решает такие уравнения, не вычисляя инверсию матрицы.

Несмотря на то, что это не стандартное математическое обозначение, MATLAB использует терминологию деления, знакомую в скалярном случае, чтобы описать решение общей системы одновременных уравнений. Два символа деления, наклонная черта,/, и наклонная черта влево, \, соответствуют этим двум функциям MATLAB mrdivide и mldivide. Эти операторы используются для двух ситуаций, где неизвестная матрица появляется слева или право на матрицу коэффициентов:

Думайте о “делении” обеих сторон уравнения Ax = b или xA = b A. Матрица коэффициентов A всегда находится в “знаменателе”.

Условия совместимости размерности для x = A\b требуют, чтобы эти две матрицы A и b имели одинаковое число строк. Решение x затем имеет одинаковое число столбцов как b и его размерность строки, равно размерности столбца A. Для x = b/A обмениваются ролями строк и столбцов.

На практике линейные уравнения формы Ax = b происходят более часто, чем те из формы xA = b. Следовательно, наклонная черта влево используется намного более часто, чем наклонная черта. Остаток от этого раздела концентрируется на операторе наклонной черты влево; соответствующие свойства оператора наклонной черты могут быть выведены из идентичности:

(b/A)' = (A'\b').

Матрица коэффициентов A не должна быть квадратной. Если A имеет размер m-by-n, то существует три случая:

Общее решение

Общее решение системы линейных уравнений Ax = b описывает все возможные решения. Можно найти общее решение:

Можно затем записать любое решение Ax = b как сумма конкретного решения Ax =b, от шага 2, плюс линейная комбинация базисных векторов от шага 1.

Остальная часть этого раздела описывает, как использовать MATLAB, чтобы найти конкретное решение Ax =b, как на шаге 2.

Квадратные системы

Наиболее распространенная ситуация включает квадратную матрицу коэффициентов A и один правый вектор-столбец стороны b.

A = pascal(3);
u = [3; 1; 4];
x = A\u

x =
      10
     -12
       5

b = magic(3);
X = A\b

X =
      19    -3    -1
     -17     4    13
       6     0    -6

P = pinv(A)*b

A = [ 1     3     7
     -1     4     4
      1    10    18 ]

rank(A)

ans =

     2

pinv(A)*b

ans =
    0.3850
   -0.1103
    0.7066

A*pinv(A)*b

ans =
    5.0000
    2.0000
   12.0000

A*pinv(A)*b

ans =
   -1.0000
    4.0000
    2.0000

rref([A b])
ans =
    1.0000         0    2.2857         0
         0    1.0000    1.5714         0
         0         0         0    1.0000

Сверхрешительные системы

Этот пример показывает, как со сверхрешительными системами часто сталкиваются в различных видах аппроксимирования кривыми к экспериментальным данным.

Количество, y, измеряется в нескольких различных значениях времени, t, чтобы произвести следующие наблюдения. Можно ввести данные и просмотреть их в таблице со следующими операторами.

t = [0 .3 .8 1.1 1.6 2.3]';
y = [.82 .72 .63 .60 .55 .50]';
B = table(t,y)

B=6×2 table
     t      y  
    ___    ____

      0    0.82
    0.3    0.72
    0.8    0.63
    1.1     0.6
    1.6    0.55
    2.3     0.5

Попытайтесь моделировать данные с затухающей показательной функцией

$y(t)=c_1+c_2{e}^{-t}$.

Предыдущее уравнение говорит, что векторный y должен быть аппроксимирован линейной комбинацией двух других векторов. Каждый - постоянный вектор, содержащий все единицы, и другой вектор с компонентами exp(-t). Неизвестные коэффициенты, $$c_1$$ и $$c_2$$, может быть вычислен путем выполнения подгонки наименьших квадратов, которая минимизирует сумму квадратов отклонений данных из модели. Существует шесть уравнений в двух неизвестных, представленных 6 2 матрица.

E = [ones(size(t)) exp(-t)]

E = 6×2

    1.0000    1.0000
    1.0000    0.7408
    1.0000    0.4493
    1.0000    0.3329
    1.0000    0.2019
    1.0000    0.1003

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

c = E\y

c = 2×1

    0.4760
    0.3413

Другими словами, подгонка наименьших квадратов к данным

Следующие операторы оценивают модель в расположенном с равными интервалами шаге в t, и затем строят результат вместе с исходными данными:

T = (0:0.1:2.5)';
Y = [ones(size(T)) exp(-T)]*c;
plot(T,Y,'-',t,y,'o')

E*c не точно равен y, но различием могут хорошо быть меньше, чем погрешности измерения в исходных данных.

Прямоугольный матричный A имеет неполный ранг, если он не имеет линейно независимых столбцов. Если A имеет неполный ранг, то решение методом наименьших квадратов к AX = B не уникально. A\B выдает предупреждение, если A имеет неполный ранг и производит решение методом наименьших квадратов. Можно использовать lsqminnorm, чтобы найти решение X, который имеет минимальную норму среди всех решений.

Недоопределенные системы

Этот пример показывает, как решение недоопределенных систем не уникально. Недоопределенные линейные системы включают больше неизвестных, чем уравнения. Матрица уехала, операция деления в MATLAB находит основное решение методом наименьших квадратов, которое имеет в большей части m ненулевые компоненты для m-by-n матрица коэффициентов.

Вот небольшой, случайный пример:

R = [6 8 7 3; 3 5 4 1]
rng(0);
b = randi(8,2,1)

R =

       6              8              7              3       
       3              5              4              1       


b =

       7       
       8      

Линейная система Rp = b вовлекает два уравнения в четыре неизвестные. Поскольку матрица коэффициентов содержит маленькие целые числа, уместно использовать команду format, чтобы отобразить решение в рациональном формате. Конкретное решение получено с

format rat
p = R\b

p =

       0       
      17/7     
       0       
     -29/7    

Одним из ненулевых компонентов является p(2), потому что R(:,2) является столбцом R с самой большой нормой. Другим ненулевым компонентом является p(4), потому что R(:,4) доминирует после того, как R(:,2) устраняется.

Полное общее решение недоопределенной системы может быть охарактеризовано путем добавления p в произвольную линейную комбинацию пустых векторов пробела, которые могут быть найдены с помощью функции null с опцией, запрашивающей рациональное основание.

Z = null(R,'r')

Z =

      -1/2           -7/6     
      -1/2            1/2     
       1              0       
       0              1       

Можно подтвердить, что R*Z является нулем и что остаточный R*x - b является маленьким для любого векторного x, где

x = p + Z*q

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

Чтобы проиллюстрировать, выберите произвольный q и создайте x.

q = [-2; 1];
x = p + Z*q;

Вычислите норму невязки.

format short
norm(R*x - b)

ans =

   2.6645e-15

Когда бесконечно много решений доступны, решение с минимальной нормой особенно интересно. Можно использовать lsqminnorm, чтобы вычислить решение методом наименьших квадратов минимальной нормы. Это решение имеет наименьшее значение для norm(p).

p = lsqminnorm(R,b)

p =

    -207/137   
     365/137   
      79/137   
    -424/137  

Решение для нескольких правых сторон

Некоторые проблемы касаются решения линейных систем, которые имеют ту же матрицу коэффициентов A, но различные правые стороны b. Когда различные значения b доступны одновременно, можно создать b как матрицу с несколькими столбцами и решить все системы уравнений одновременно с помощью одной команды наклонной черты влево: X = A\[b1 b2 b3 …].

Однако иногда различные значения b не все доступны одновременно, что означает, что необходимо решить несколько систем уравнений последовательно. Когда вы решаете одну из этих систем уравнений с помощью наклонной черты (/) или наклонная черта влево (\), оператор разлагает на множители матрицу коэффициентов A и использует это матричное разложение, чтобы вычислить решение. Однако в каждый последующий раз, когда вы решаете аналогичную систему уравнений с различным b, оператор вычисляет то же разложение A, который является избыточным вычислением.

Решение этой проблемы состоит в том, чтобы предварительно вычислить разложение A, и затем снова использовать факторы, чтобы решить для различных значений b. На практике, однако, предварительное вычисление разложения этим способом может быть трудным, поскольку необходимо знать, какое разложение вычислить (LU, LDL, Холесский, и так далее), а также как умножить факторы, чтобы решить проблему. Например, с разложением LU необходимо решить две линейных системы, чтобы решить исходную систему Ax = b:

[L,U] = lu(A);
x = U \ (L \ b);

Вместо этого рекомендуемый метод для решения линейных систем с несколькими последовательными правыми сторонами должен использовать объекты decomposition. Эти объекты позволяют вам усилить выигрыши в производительности предварительного вычисления матричного разложения, но они не требуют знания того, как использовать матричные факторы. Можно заменить предыдущее разложение LU на:

dA = decomposition(A,'lu');
x = dA\b;

Если вы не уверены, какое разложение использовать, decomposition(A) выбирает правильный тип на основе свойств A, подобного тому, что делает наклонная черта влево.

Вот простой тест возможных выигрышей в производительности этого подхода. Тест решает ту же разреженную линейную систему 100 раз с помощью и наклонной черты влево (\) и decomposition.

n = 1e3;
A = sprand(n,n,0.2) + speye(n);
b = ones(n,1);

% Backslash solution
tic
for k = 1:100
    x = A\b;
end
toc

Elapsed time is 9.006156 seconds.

% decomposition solution
tic
dA = decomposition(A);
for k = 1:100
    x = dA\b;
end
toc

Elapsed time is 0.374347 seconds.

Для этой проблемы решение decomposition намного быстрее, чем использование одной только наклонной черты влево, все же синтаксис остается простым.

Итеративные методы

Если матрица коэффициентов A является большой и разреженной, методы факторизации обычно не эффективны. Iterative methods генерирует серию приближенных решений. MATLAB обеспечивает несколько итерационных методов, чтобы обработать большие, разреженные входные матрицы.

Многопоточное вычисление

MATLAB поддерживает многопоточное вычисление во многой линейной алгебре и поэлементных числовых функциях. Эти функции автоматически выполняются на нескольких потоках. Для функции или выражения, чтобы выполниться быстрее на нескольких центральных процессорах, много условий должны быть верными:

inv, lscov, linsolve и mldivide показывают значительное увеличение скорости на больших массивах с двойной точностью (на порядке 10 000 элементов или больше), когда многопоточность включена.