LDLT-разложение

[nabbla1]

В некотором смысле, "заметки для себя": алгоритмы линейной алгебры обладают нехорошим свойством write-only, особенно если там участвует "треугольная адресация". Пока разбираешь - всё вроде понятно, а чуть-чуть погодя смотришь на это безобразие - и оно превращается в набор вложенных циклов и десяток индексов массива, которые обновляются на каждой

Comments

[lj_frank_bot]

Здравствуйте!
Система категоризации Живого Журнала посчитала, что вашу запись можно отнести к категории: Работа.
Если вы считаете, что система ошиблась - напишите об этом в ответе на этот комментарий. Ваша обратная связь поможет сделать систему точнее.
Фрэнк,
команда ЖЖ.

[nabbla1]

Наука.

[lj_frank_bot]

Упс.

[ddem1979]

Позанудствую, но на телефоне андроидный файрфокс на половину формул выдал кашу...
\left он не понимает...

[nabbla1]

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

Вижу две возможные причины: или была ошибка на серверной стороне (сайт latex.codecogs.com), или выдаваемый .svg браузер невзлюбил.

Там можно в .png результат выдавать, но мне .svg больше нравился, т.к при масштабировании странички формулы остаются гладкими. Давайте сейчас переправлю в .png - и посмотрите ещё разок. (там просто вместо latex.codecogs.com/svg.latex?... пишется latex.codecogs.com/png.latex?...)

[ddem1979]

Вот оно сейчас:
https://cloud.mail.ru/public/2Jgv/3HgWP4kig

[ddem1979]

Я ТеХ люблю, но и проблем с ним много. Нормальный шрифт и пробелы в формулах. оформить все по ГОСТам, которые каждый год меняются. В-общем, я в ТеХе только научное набираю, на отсылку наружу, а всю повседневку - в Ворд 2003 (((

[kouzdra]

алгоритмы линейной алгебры обладают нехорошим свойством write-only, особенно если там участвует "треугольная адресация". Пока разбираешь - всё вроде понятно, а чуть-чуть погодя смотришь на это безобразие - и оно превращается в набор вложенных циклов и десяток индексов массива, которые обновляются на каждой итерации не вполне понятным образом

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

[nabbla1]

Я пока только встречал библиотеку на C, которая была сделана честной переделкой кода на фортране, с индексами типа ii, jj, kk, nn, mm. Возможно, решили "работает-не лезь!".

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

[kouzdra]

Я скорее про что-то APL-образное, только подточенное именно под матрицы-вектора.

[quartisz]

А где можно найти готовую процедуру для LDLt на Pascal или Фортран, ну или на Си на худой конец?

[nabbla1]

У меня есть:

procedure CholeskyLDL (const a: array of Real; N: Integer; out u: array of Real);
var i: Integer; //счётчик для самого глубокого цикла
icol: Integer; //столбец, на котором текущий обраб. элемент
irow: Integer; //строка на которой текущий обраб. элемент
j: Integer; //указывает на элемент u[1, icol]
k: Integer; //указывает на элемент u[irow, icol], который и обрабатывается
l: Integer; //указывает на элемент u[i, irow]
m: Integer; //указывает на элемент u[i, icol]
d: Integer; //указывает на диагональный элемент u[i,i]
w: Real; //промежут. значение

begin
if N <= 0 then
raise Exception.Create('CholeskyLDL: N <= 0');

if ( Length(a) < N * ( N + 1 ) div 2 ) or (Length(u) < N * (N + 1) div 2) then

[quartisz]

Спасибо. А массивы почему одномерные?
Я вот такую написал процедуру.

type
TExtVec = array of Extended;
TExtMatr = array of TExtVec;

procedure LDLT(n: Integer; var a,L: TExtMatr; var d: TExtVec);
var i,j,k: Integer;
s: Extended;
begin
SetLength(L,n+1,n+1);
SetLength(d,n+1);
s:= 0;
for i:= 1 to n do
for j:= i+1 to n do
begin
s:= A[j][i];
for k:= 1 to i do
s:= s-L[i,k]*D[k]* L[j,k];
if (i = j) then
begin
D[i]:= s;//диагональный элемент
L[i,i]:= 1;//диагональ
end
else L[j,i]:= s/D[i];//внедиагональный элемент
end;
end;

[nabbla1]

Я брал пример с этого дела: https://people.sc.fsu.edu/~jburkardt/f77_src/asa007/asa007.html

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

[ext_6383178]

Вопрос от дилетанта: разложение Холецкого на множители LDLt позволяет избежать вычисления квадратных корней, делает ли это возможным применять данное разложение в отношении матрицы, не являющейся положительно определенной?

[nabbla1]

Квадратные корни здесь не нужны, да.

Само по себе LDLT-разложение допустимо для неотрицательно определённой матрицы, но я не проверял, сможет ли приведённый алгоритм это сделать, или "поделит на ноль" по чём зря. Возможно, придётся этот случай отдельно рассмотреть.

И пишут, что и более общего вида матрицы можно так раскладывать, главное, чтобы все её главные миноры были ненулевыми. Но и это не пробовал, не было необходимости.