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

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

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

В матричном обозначении общая проблема принимает следующую форму: Учитывая две матрицы A и b, действительно там существует уникальная матрица x, так, чтобы Ax = b или xA = 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 на n размера, то существует три случая:

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

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

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

Остальная часть этого раздела описывает, как использовать MATLAB, чтобы найти конкретное решение Оси =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 должен быть аппроксимирован линейной комбинацией двух других векторов. Каждый - постоянный вектор, содержащий все единицы, и другой вектор с компонентами exp(-t). Неизвестные коэффициенты, и, могут быть вычислены путем выполнения соответствия наименьших квадратов, которое минимизирует сумму квадратов отклонений данных от модели. Существует шесть уравнений в двух неизвестных, представленных 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 является большой и разреженной, методы факторизации обычно не эффективны. Итеративные методы генерируют серию приближенных решений. MATLAB предоставляет несколько итеративных методов, чтобы обработать большие, разреженные входные матрицы.

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

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

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