18.3 Разложение матриц

: R = chol (A)

: [R, p] = chol (A)

: [R, p, Q] = chol (A)

: [R, p, Q] = chol (A, "vector")

: [L, …] = chol (…, "lower")

: [R, …] = chol (…, "upper")

Вычислить верхний фактор Холецкого, R, вещественной симметричной или комплексной эрмитовой положительно определённой матрицы A.

Верхний фактор Холецкого R вычисляется, используя верхнюю треугольную часть матрицы A и определяется как

R' * R = A.

Вызов chol с необязательным "upper" флагом имеет тот же результат. В отличие от этого, с необязательным флагом "lower", chol возвращает нижнее треугольное разложение, вычисленное с использованием нижней треугольной части матрицы A, так что

L * L' = A.

Вызов с одним выходным аргументом chol завершается ошибкой, если матрица A не является положительно определённой. Обратите внимание, что если матрица A не является вещественной симметричной или комплексной эрмитовой, то нижняя треугольная часть рассматривается как (комплексно-сопряжённый) транспонированная верхняя треугольная часть, или наоборот, в зависимости от "lower" флага.

Вызов с двумя или более выходными аргументами p указывает, была ли матрица A положительно определённой, и chol не завершается ошибкой. Значение 0 для p означает, что матрица A является положительно определённой, и R даёт разложение. В противном случае p будет иметь положительное значение.

Если вызывается с тремя выходными аргументами, матрица A должна быть разреженной, и перед факторизацией к матрице A применяется сохраняющая разреженность перестановка строк/столбцов. То есть R является разложением A(Q,Q) таким образом

R' * R = Q' * A * Q.

Сохраняющая разреженность перестановка обычно возвращается как матрица. Однако, с необязательным флагом "vector", Q будет возвращена как вектор, так что

R' * R = A(Q, Q).

В общем случае, нижнее треугольное разложение значительно быстрее для разреженных матриц.

См. также: hess, lu, qr, qz, schur, svd, ichol, cholinv, chol2inv, cholupdate, cholinsert, choldelete, cholshift.

: cholinv (A)

Вычислить обратную симметричную положительно определённую матрицу A с использованием разложения Холецкого.

См. также: chol, chol2inv, inv.

: chol2inv (U)

Найти обратную симметричной положительно определённой квадратной матрице по её разложению Холецкого, U.

Обратите внимание, что U должна быть верхней треугольной матрицей с положительными диагональными элементами. chol2inv (U) предоставляет inv (U'*U), но это намного быстрее, чем использование inv.

См. также: chol, cholinv, inv.

: [R1, info] = cholupdate (R, u, op)

Обновление или понижение разложения Холецкого.

При заданной верхней треугольной матрице R и столбце-векторе u, попытаться определить другую верхнюю треугольную матрицу R1, такую что

• R1’*R1 = R’*R + u*u’ если op — "+"
• R1’*R1 = R’*R - u*u’ если op — "-"

Если op — "-", info устанавливается в

• 0, если понижение прошло успешно,
• 1, если R’*R - u*u’ не является положительно определённой,
• 2, если R вырождена.

Если info не указан, сообщение об ошибке выводится в случаях 1 и 2.

См. также: chol, cholinsert, choldelete, cholshift.

: R1 = cholinsert (R, j, u)

: [R1, info] = cholinsert (R, j, u)

Обновление разложения Холецкого при вставке строки или столбца в исходную факторизованную матрицу.

При заданном разложении Холецкого вещественной симметричной или комплексной эрмитовой положительно определённой матрицы A = R’*R, где R верхняя треугольная, вернуть разложение Холецкого A1, где A1(p,p) = A, A1(:,j) = A1(j,:)’ = u и p = [1:j-1,j+1:n+1]. u(j) должно быть положительным.

При возвращении, info устанавливается в

• 0, если вставка прошла успешно,
• 1, если A1 не является положительно определённой,
• 2, если R вырождена.

Если info не указан, сообщение об ошибке выводится в случаях 1 и 2.

См. также: chol, cholupdate, choldelete, cholshift.

: R1 = choldelete (R, j)

Обновление разложения Холецкого с удалением строки или столбца из исходной факторизованной матрицы.

При заданном разложении Холецкого вещественной симметричной или комплексной эрмитовой положительно определённой матрицы A = R’*R, где R верхняя треугольная, вернуть разложение Холецкого A(p,p), где p = [1:j-1,j+1:n+1].

См. также: chol, cholupdate, cholinsert, cholshift.

: R1 = cholshift (R, i, j)

Обновление разложения Холецкого с изменением диапазона столбцов в исходной факторизованной матрице.

При заданном разложении Холецкого вещественной симметричной или комплексной эрмитовой положительно определённой матрицы A = R’*R, где R верхняя треугольная, вернуть разложение Холецкого A(p,p), где p — перестановка
p = [1:i-1, shift(i:j, 1), j+1:n] если i < j
или
p = [1:j-1, shift(j:i,-1), i+1:n] если j < i.

См. также: chol, cholupdate, cholinsert, choldelete.

: H = hess (A)

: [P, H] = hess (A)

Вычислить разложение Хессенберга матрицы A.

Разложение Хессенберга — это P * H * P' = A где P — квадратная унитарная матрица (P' * P = I, используя комплексно-сопряжённое транспонирование) и H — верхняя матрица Хессенберга (H(i, j) = 0 forall i > j+1).

Разложение Хессенберга обычно используется как первый шаг в вычислении собственных значений, но также имеет и другие применения (см. Golub, Nash, and Van Loan, IEEE Transactions on Automatic Control, 1979).

См. также: eig, chol, lu, qr, qz, schur, svd.

: [L, U] = lu (A)

: [L, U, P] = lu (A)

: [L, U, P, Q] = lu (S)

: [L, U, P, Q, R] = lu (S)

: […] = lu (S, thresh)

: y = lu (…)

: […] = lu (…, "vector")

Вычислить LU-разложение матрицы A.

Если A — полная матрица, то используются подпрограммы из LAPACK, а если A — разреженная, то используется UMFPACK.

Результат возвращается в переупорядоченном виде в соответствии с необязательным возвращаемым значением P. Например, для матрицы A = [1, 2; 3, 4],

[L, U, P] = lu (A)

возвращается

L =

  1.00000  0.00000
  0.33333  1.00000

U =

  3.00000  4.00000
  0.00000  0.66667

P =

  0  1
  1  0

Матрица не обязана быть квадратной.

При вызове с двумя или тремя выходными аргументами и разреженной входной матрицей, lu не пытается выполнить сохраняющие разреженность перестановки столбцов. При вызове с четвёртым выходным аргументом, возвращается сохраняющая разреженность преобразование столбцов Q, такое что P * A * Q = L * U. Это предпочтительный способ вызова lu с разреженными входными матрицами.

При вызове с пятым выходным аргументом и разреженной входной матрицей, lu пытается использовать множитель масштабирования R для входной матрицы таким образом, что P * (R \ A) * Q = L * U. Это обычно приводит к более разреженному и более стабильному разложению.

Можно задать дополнительный входной аргумент thresh, определяющий порог для выбора опорного элемента. thresh может быть скаляром, в этом случае он определяет порог выбора опорного элемента в UMFPACK для симметричных и несимметричных случаев. Если thresh является 2-элементным вектором, то первый элемент определяет порог выбора опорного элемента для несимметричной стратегии выбора опорного элемента в UMFPACK, а второй — для симметричной стратегии. По умолчанию используются значения, определённые spparms ([0.1, 0.001]).

При задании строкового аргумента "vector", lu возвращает значения P и Q как векторные значения, так что для полной матрицы A(P,:) = L * U, и R(P,:) * A(:,Q) = L * U.

С двумя выходными аргументами возвращает переупорядоченные формы верхней и нижней треугольных матриц, такие что A = L * U. С одним выходным аргументом y возвращается матрица, возвращаемая процедурами LAPACK. Если входная матрица является разреженной, то матрица L встраивается в U для получения значения, аналогичного полному случаю. Для полных и разреженных матриц, lu теряет информацию о перестановке.

См. также: luupdate, ilu, chol, hess, qr, qz, schur, svd.

: [L, U] = luupdate (L, U, x, y)

: [L, U, P] = luupdate (L, U, P, x, y)

Дано LU-разложение вещественной или комплексной матрицы A = L*U, где L — нижняя единичная трапециевидная матрица, а U — верхняя трапециевидная матрица. Возвращает LU-разложение A + x*y’, где x и y — столбцовые векторы (обновление ранга 1) или матрицы с одинаковым числом столбцов (обновление ранга k).

Необязательно, можно использовать обновление с выбором строк, предоставив матрицу перестановок строк (P); в этом случае возвращается обновлённая матрица перестановок. Обратите внимание, что если L, U, P — LU-разложение с выбором строк, полученное с помощью lu:

[L, U, P] = lu (A);

то разложение A+x*y.' можно получить либо как

[L1, U1] = lu (L, U, P*x, y)

либо как

[L1, U1, P1] = lu (L, U, P, x, y)

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

Случай матрицы выполняется как последовательность обновлений ранга 1; следовательно, для достаточно большого k будет и быстрее, и точнее пересчитать разложение с нуля.

См. также: lu, cholupdate, qrupdate.

: [Q, R] = qr (A)

: [Q, R, P] = qr (A) # non-sparse A

: X = qr (A) # non-sparse A

: R = qr (A) # sparse A

: [C, R] = qr (A, B)

: […] = qr (…, 0)

: […] = qr (…, "vector")

: […] = qr (…, "matrix")

Вычислить QR-разложение A, используя стандартные подпрограммы LAPACK.

QR-разложение имеет вид

Q * R = A

где Q – ортогональная матрица, а R – верхняя треугольная матрица.

Например, для матрицы A = [1, 2; 3, 4],

[Q, R] = qr (A)

возвращает

Q =

  -0.31623  -0.94868
  -0.94868   0.31623

R =

  -3.16228  -4.42719
   0.00000  -0.63246

что при перемножении возвращает исходную матрицу

Q * R
  ⇒
     1.0000   2.0000
     3.0000   4.0000

Если запрашивается только одно значение, то это либо R, если A разреженная, либо X, такое что R = triu (X) , если A полная. (Примечание: в отличие от большинства команд, единственное возвращаемое значение не является первым, когда запрашивается несколько значений.)

Если матрица A полная, и запрашивается третий выход P, то qr вычисляет переупорядоченное QR-разложение

Q * R = A * P

где Q – ортогональная матрица, R – верхняя треугольная матрица, а P – матрица перестановок.

Переупорядоченное QR-разложение обладает дополнительным свойством, что диагональные элементы R упорядочены по убыванию величины. Другими словами, abs (diag (R)) будут отсортированы от наибольшего к наименьшему.

Например, для матрицы A = [1, 2; 3, 4],

[Q, R, P] = qr (A)

возвращает

Q =

  -0.44721  -0.89443
  -0.89443   0.44721

R =

  -4.47214  -3.13050
   0.00000   0.44721

P =

   0  1
   1  0

Если входная матрица A разреженная, то разреженное QR-разложение вычисляется с помощью CSPARSE. Поскольку матрица Q в общем случае является полной матрицей, рекомендуется запрашивать только одно возвращаемое значение R. В этом случае вычисление избегает построения Q и возвращает R такое, что R = chol (A' * A).

Если дополнительно задана матрица B и запрашиваются два возвращаемых значения, то qr возвращает C, где C = Q' * B. Это позволяет вычислить приближение наименьших квадратов A \ B как

[C, R] = qr (A, B)
x = R \ C

Если последним аргументом является строка "vector", то P является вектором перестановок (столбцов A), а не матрицей перестановок. В этом случае определяющее соотношение имеет вид

Q * R = A(:, P)

Однако по умолчанию возвращается матрица перестановок, и это можно явно указать, используя в качестве последнего аргумента "matrix".

Если последним аргументом является скаляр 0, то возвращается "economy" разложение. Когда исходная матрица A имеет размер MxN и M > N, то "economy" разложение вычислит только N строк в R и N столбцов в Q и опуская нули в R. Если M ≤ N, то нет разницы между экономическим и стандартным разложениями. При вычислении "economy" разложения выход P всегда является вектором, а не матрицей.

Предыстория: QR-разложение применяется при решении задач наименьших квадратов

min norm (A*x - b)

для переопределённых систем уравнений (т.е., A – высокая, узкая матрица).

Переупорядоченное QR-разложение [Q, R, P] = qr (A) позволяет построить ортогональную базу span (A).

См. также: chol, hess, lu, qz, schur, svd, qrupdate, qrinsert, qrdelete, qrshift.

: [Q1, R1] = qrupdate (Q, R, u, v)

Обновить QR-разложение, используя векторы или матрицы обновления.

Дано QR-разложение вещественной или комплексной матрицы A = Q*R, где Q унитарная, а R верхняя трапециевидная, вернуть QR-разложение A + u*v’, где u и v – столбцовые векторы (обновление ранга 1) или матрицы с одинаковым числом столбцов (обновление ранга k). Обратите внимание, что последний случай выполняется как последовательность обновлений ранга 1; таким образом, при достаточно большом k будет быстрее и точнее перевычислить разложение с нуля.

Поставленное QR-разложение может быть полным (Q квадратная) или экономичным (R квадратная).

См. также: qr, qrinsert, qrdelete, qrshift.

: [Q1, R1] = qrinsert (Q, R, j, x, orient)

Обновление QR-разложения, учитывая строку или столбец для вставки в исходную матрицу.

Дано QR-разложение вещественной или комплексной матрицы A = Q*R, где Q унитарна, а R верхняя трапециевидная. Возвращает QR-разложение матрицы [A(:,1:j-1) x A(:,j:n)], где u — столбец, который нужно вставить в A (если orient равно "col"), или QR-разложение матрицы [A(1:j-1,:);x;A(:,j:n)], где x — строка, которую нужно вставить в A (если orient равно "row").

Значение по умолчанию для orient равно "col". Если orient равно "col", u может быть матрицей, а j — вектором индексов, что приводит к QR-разложению матрицы B, такой что B(:,j) даёт u, а B(:,j) = [] даёт A. Обратите внимание, что последний случай выполняется как последовательность k вставок; таким образом, при достаточно большом k будет быстрее и точнее повторно вычислить разложение с нуля.

Если orient равно "col", предоставленное QR-разложение может быть полным (Q квадратная) или экономичным (R квадратная).

Если orient равно "row", необходимо полное разложение.

См. также: qr, qrupdate, qrdelete, qrshift.

: [Q1, R1] = qrdelete (Q, R, j, orient)

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

Дано QR-разложение вещественной или комплексной матрицы A = Q*R, где Q унитарна, а R верхняя трапециевидная. Возвращает QR-разложение матрицы [A(:,1:j-1), U, A(:,j:n)], где u — столбец, который нужно вставить в A (если orient равно "col"), или QR-разложение матрицы [A(1:j-1,:);X;A(:,j:n)], где x — строка (orient равно "row"). Значение по умолчанию для orient равно "col".

Если orient равно "col", j может быть вектором индексов, что приводит к QR-разложению матрицы B, такой что A(:,j) = [] даёт B. Обратите внимание, что последний случай выполняется как последовательность k удалений; таким образом, при достаточно большом k будет быстрее и точнее повторно вычислить разложение с нуля.

Если orient равно "col", предоставленное QR-разложение может быть полным (Q квадратная) или экономичным (R квадратная).

Если orient равно "row", необходимо полное разложение.

См. также: qr, qrupdate, qrinsert, qrshift.

: [Q1, R1] = qrshift (Q, R, i, j)

Обновление QR-разложения, учитывая диапазон столбцов для сдвига в исходной матрице.

Дано QR-разложение вещественной или комплексной матрицы A = Q*R, где Q унитарна, а R верхняя трапециевидная. Возвращает QR-разложение матрицы A(:,p), где p — перестановка
p = [1:i-1, shift(i:j, 1), j+1:n] если i < j
или
p = [1:j-1, shift(j:i,-1), i+1:n] если j < i.

См. также: qr, qrupdate, qrinsert, qrdelete.

: lambda = qz (A, B)

: [AA, BB, Q, Z, V, W, lambda] = qz (A, B)

: [AA, BB, Z] = qz (A, B, opt)

: [AA, BB, Z, lambda] = qz (A, B, opt)

Вычисление QZ-разложения задачи обобщённых собственных значений.

Задача обобщённых собственных значений определяется как

A x = lambda B x

Существует три формы вызова функции:

1. lambda = qz (A, B)

   Вычисление обобщённых собственных значений lambda.

2. [AA, BB, Q, Z, V, W, lambda] = qz (A, B)

   Вычисление QZ-разложения, обобщённых собственных векторов и обобщённых собственных значений.

   AA = Q * A * Z, BB = Q * B * Z
   A * V = B * V * diag (lambda)
   W' * A = diag (lambda) * W' * B

   с Q и Z ортогональными (унитарными для комплексного случая).

3. [AA, BB, Z {, lambda}] = qz (A, B, opt)

   Как в форме 2 выше, но позволяет упорядочить обобщённые собственные пары для, например, решения дискретных алгебраических уравнений Риккати. Форма 3 недоступна для комплексных матриц и не вычисляет обобщённые собственные векторы V, W, ни ортогональную матрицу Q.

   opt

   для упорядочения собственных значений обобщённой задачи собственных значений. Ведущий блок пересмотренного разложения содержит все собственные значения, которые удовлетворяют:

   "N"

   неупорядоченные (по умолчанию)

   "S"

   малые: ведущий блок содержит все |lambda| < 1

   "B"

   большие: ведущий блок содержит все |lambda| ≥ 1

   "-"

   отрицательная вещественная часть: ведущий блок содержит все собственные значения в открытой левой полуплоскости

   "+"

   неотрицательная вещественная часть: ведущий блок содержит все собственные значения в замкнутой правой полуплоскости

Примечание: qz выполняет перестановку балансировки, но не масштабирование (см. balance), что может привести к менее точным результатам, чем eig. Порядок выходных аргументов был выбран для совместимости с MATLAB.

См. также: eig, gsvd, balance, chol, hess, lu, qr, qzhess, schur.

: [aa, bb, q, z] = qzhess (A, B)

Вычисление разложения Гессенберга-треугольника матричного разложения (A, B), возвращающее aa = q * A * z, bb = q * B * z, с q и z ортогональными.

Например:

[aa, bb, q, z] = qzhess ([1, 2; 3, 4], [5, 6; 7, 8])
  ⇒ aa =
      -3.02244  -4.41741
       0.92998   0.69749
  ⇒ bb =
      -8.60233  -9.99730
       0.00000  -0.23250
  ⇒ q =
      -0.58124  -0.81373
      -0.81373   0.58124
  ⇒ z =
     Diagonal Matrix
       1   0
       0   1

Разложение Гессенберга-треугольника — это первый шаг в алгоритме QZ-разложения Молера и Стюарта.

Алгоритм взят из Голуба и Ван Лоана, Вычисления с матрицами, 2-е издание.

См. также: lu, chol, hess, qr, qz, schur, svd.

: S = schur (A)

: S = schur (A, "real")

: S = schur (A, "complex")

: S = schur (A, opt)

: [U, S] = schur (…)

Вычисление разложения Шюра матрицы A.

Разложение Шюра определяется как

S = U' * A * U

где U — унитарная матрица (U'* U — единичная матрица) и S — верхняя треугольная матрица. Собственные значения A (и S) — это диагональные элементы S. Если матрица A вещественная, то вычисляется вещественное разложение Шюра, в котором матрица U ортогональна, а S — блочно-верхняя треугольная с блоками размера не более 2 x 2 по диагонали. Диагональные элементы S (или собственные значения 2 x 2 блоков, при необходимости) — собственные значения A и S.

По умолчанию для вещественных матриц используется вещественное разложение Шюра. Комплексное разложение можно принудительно задать, передав флаг "complex".

Собственные значения необязательно упорядочены по диагонали в соответствии со значением opt. opt = "a" указывает, что все собственные значения с отрицательной вещественной частью должны быть перемещены в ведущий блок S (используется в are), opt = "d" указывает, что все собственные значения с модулем меньше единицы должны быть перемещены в ведущий блок S (используется в dare), а opt = "u", по умолчанию, указывает, что никакого упорядочения собственных значений не должно происходить. Ведущие k столбцов U всегда охватывают A-инвариантное подпространство, соответствующее k ведущим собственным значениям S.

Разложение Шюра используется для вычисления собственных значений квадратной матрицы и имеет приложения в решении алгебраических уравнений Риккати в управлении (см. are и dare).

См. также: rsf2csf, ordschur, ordeig, lu, chol, hess, qr, qz, svd.

: [U, T] = rsf2csf (UR, TR)

Преобразование вещественной, верхней квазитреугольной формы Шюра TR в комплексную, верхнюю треугольную форму Шюра T.

Обратите внимание, что следующие соотношения выполняются:

UR * TR * UR' = U * T * U' и U' * U — единичная матрица I.

Также обратите внимание, что U и T не уникальны.

См. также: schur.

: [UR, SR] = ordschur (U, S, select)

Переупорядочивает факторизацию по реальному Шуру (U,S), полученную с помощью функции schur, так, чтобы выбранные собственные значения появлялись в диагональных блоках левого верхнего угла квазитреугольной матрицы Шура.

Логический вектор select определяет выбранные собственные значения, так как они появляются вдоль диагонали S.

Например, задана матрица A = [1, 2; 3, 4], и ее разложение Шура

[U, S] = schur (A)

которое возвращает

U =

  -0.82456  -0.56577
   0.56577  -0.82456

S =

  -0.37228  -1.00000
   0.00000   5.37228

Возможно переупорядочить разложение так, чтобы положительное собственное значение было в левом верхнем углу, выполнив:

[U, S] = ordschur (U, S, [0,1])

См. также: schur, ordeig.

: lambda = ordeig (A)

: lambda = ordeig (A, B)

Возвращает собственные значения квазитреугольных матриц в порядке их появления в матрице A.

Квазитреугольная матрица A обычно является результатом факторизации Шура. Если вызвана со вторым аргументом B, то обобщенные собственные значения пары A, B возвращаются в порядке появления матрицы A-lambda*B. Пара A, B обычно является результатом разложения QZ.

См. также: ordschur, eig, schur, qz.

: angle = subspace (A, B)

Определяет наибольший главный угол между двумя подпространствами, охватываемыми столбцами матриц A и B.

: s = svd (A)

: [U, S, V] = svd (A)

: [U, S, V] = svd (A, "econ")

: [U, S, V] = svd (A, 0)

Вычислить сингулярное разложение A.

Сингулярное разложение определяется соотношением

A = U*S*V'

Функция svd обычно возвращает только вектор сингулярных значений. При вызове с тремя значениями возврата она вычисляет U, S и V. Например,

svd (hilb (3))

возвращает

ans =

  1.4083189
  0.1223271
  0.0026873

и

[u, s, v] = svd (hilb (3))

возвращает

u =

  -0.82704   0.54745   0.12766
  -0.45986  -0.52829  -0.71375
  -0.32330  -0.64901   0.68867

s =

  1.40832  0.00000  0.00000
  0.00000  0.12233  0.00000
  0.00000  0.00000  0.00269

v =

  -0.82704   0.54745   0.12766
  -0.45986  -0.52829  -0.71375
  -0.32330  -0.64901   0.68867

При указании второго аргумента, не равного 0, svd возвращает разложение экономичного размера, устраняя ненужные строки или столбцы U или V.

Если второй аргумент равен ровно 0, то выбор разложения основан на матрице A. Если у A больше строк, чем столбцов, возвращается разложение экономичного размера, в противном случае вычисляется обычное разложение.

Примечания по алгоритму: При вычислении полного разложения (левые и правые сингулярные матрицы в дополнение к сингулярным значениям) есть выбор из двух процедур в LAPACK. По умолчанию в Octave используется gesvd. Альтернатива - gesdd, которая в 5 раз быстрее, но может использовать больше памяти и быть неточной для некоторых входных матриц. Обратитесь к документации по svd_driver для получения дополнительной информации о выборе драйвера.

См. также: svd_driver, svds, eig, lu, chol, hess, qr, qz.

: val = svd_driver ()

: old_val = svd_driver (new_val)

: svd_driver (new_val, "local")

Запрос или установка базового драйвера LAPACK, используемого функцией svd.

В настоящее время распознаются значения "gesdd" и "gesvd". По умолчанию "gesvd".

При вызове изнутри функции с опцией "local", переменная изменяется локально для функции и любых вызываемых ею подпрограмм. Исходное значение переменной восстанавливается при выходе из функции.

Примечания по алгоритму: Библиотека LAPACK предоставляет две процедуры для вычисления полного сингулярного разложения (левые и правые сингулярные матрицы, а также сингулярные значения). При вычислении только сингулярных значений данное обсуждение не имеет отношения.

Новая процедура gesdd основана на алгоритме "Разделяй и властвуй", который в 5 раз быстрее, чем альтернативная gesvd, основанная на факторизации QR. Однако новый алгоритм может использовать значительно больше памяти. Для матрицы MxN размер памяти составляет порядка O(min(M,N) ^ 2), в то время как альтернатива - порядка O(max(M,N)).

Помимо вопросов скорости и памяти, были случаи, когда некоторые входные матрицы не были точно разложены с помощью gesdd. См. текущую открытую ошибку https://savannah.gnu.org/bugs/?55564. До решения этих проблем с точностью в новой версии библиотеки LAPACK, драйвер по умолчанию в Octave установлен на "gesvd".

См. также: svd.

: [housv, beta, zer] = housh (x, j, z)

Вычислить вектор отражения Хаусхолдера housv, чтобы отразить x, чтобы он стал j-ым столбцом единичной матрицы, т.е.

(I - beta*housv*housv')x =  norm (x)*e(j) if x(j) < 0,
(I - beta*housv*housv')x = -norm (x)*e(j) if x(j) >= 0

Входы

x

вектор

j

индекс в векторе

z

порог для нуля (обычно должен быть числом 0)

Выходы (см. Голуб и Ван Лоан):

beta

Если beta = 0, то отражение не нужно применять (zer устанавливается в 0)

housv

вектор Хаусхолдера

: [u, h, nu] = krylov (A, V, k, eps1, pflg)

Построить ортогональную базу u блока подпространства Кристоффеля.

Блок подпространства Кристоффеля имеет следующий вид:

[v a*v a^2*v … a^(k+1)*v]

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

Если V - вектор, то h содержит матрицу Хессенберга, такую что a*u == u*h+rk*ek', где rk = a*u(:,k)-u*h(:,k), а ek' - вектор [0, 0, …, 1] длины k. В противном случае h бессмысленно.

Если V - вектор и k больше length (A) - 1, то h содержит матрицу Хессенберга, такую что a*u == u*h.

Значение nu - это размерность размаха подпространства Кристоффеля (основанного на eps1).

Если b - вектор, и k больше m-1, то h содержит матрицу Хессенберга разложения A.

Необязательный параметр eps1 - порог для нуля. Значение по умолчанию равно 1e-12.

Если необязательный параметр pflg отличен от нуля, используется поворот строк для повышения численной устойчивости. Значение по умолчанию равно 0.

Ссылка: A. Hodel, P. Misra, Частичный поворот при вычислении подпространств Кристоффеля больших разреженных систем, Труды 42-й конференции IEEE по принятию решений и управлению, декабрь 2003 г.