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

Вычислительные соображения

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

В матричной нотации общая задача принимает следующий вид: Учитывая две матрицы 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-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){}_{\mathrm{=}}{}_{\mathrm{c1}}{}^{+}$c2e-t.

Предыдущее уравнение говорит, что вектор y должен быть аппроксимирован линейной комбинацией двух других векторов. Один - постоянный вектор, содержащий все, а другой - вектор с компонентами exp(-t). Неизвестные коэффициенты, $${}_{\mathrm{c1}}$$ и $${}_{\mathrm{c2}}$$, могут быть вычислены путем аппроксимации методом наименьших квадратов, что минимизирует сумму квадратов отклонений данных от модели. Существует шесть уравнений в двух неизвестных, представленных матрицей 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около-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, Cholesky и так далее), а также как умножить факторы для решения проблемы. Например, при разложении LU необходимо решить две линейные системы для решения исходной системы Ax = b:

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

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

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 решение гораздо быстрее, чем использование обратной косой черты, но синтаксис остается простым.

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

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

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

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

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