Требуется оптимизировать функцию ЦОС, в которой наиболее затратными являются функция умножения комплексных матриц и функция обратной матрицы.

Я заметил, что перемножаемые матрицы являются якобы теплицевыми, но не симметричными и не квадратными, т.е. примерно следующих форматов:

{ a+0, a+1, a+2, a+3;
a-1, a+0, a+1, a+2;
a-2, a-1, a+0, a+1 }

{ a+0, a+1, a+2;
a-1, a+0, a+1;
a-2, a-1, a+0;
a-3, a-2, a-1 }

Встречал алгоритмы, которые позволяют работать с теплицевыми матрицами размером [NxN].

Существуют ли оптимизированные алгоритмы, которые позволяют перемножать теплицевые не квадратные матрицы размером [NxM], где N != M ???
Может кто-нибудь встречал, подскажите, пожалуйста...

может быть теплицевыми являются только квадратные матрицы?

Квадратичный алгоритм описан здесь:
http://stackoverflow.com/questions/1588952...eplitz-matrices

Элементы на диагоналях получаются с помощью скользящего окна. Наверное, можно еще улучшить, считая, что умножается A * transpose(A)

Речь случайно не о методе наименьших квадратов или какой-то его разновидности? Если да, то полагаю у вас следующее: (A^H * A)^-1 * A^H. Ищется псевдоинверсия от матрицы данных А без прямого рассчета корреляционной матрицы Ф = А^H * A. В таком случае может вам лучше посмотреть в сторону svd реализации псевдоинверсии? Решение устойчивое и много где описано, например у Хайкина Adaptive filter theory. А насчет Теплица, корреляционная матрица уже не Теплицева, так что с инверсией (A^H * A)^-1 надо уже думать отдельно. Если ее оценивать по данным/матрице данных, она вроде как сходится к Эрмитовой со временем, т.к. по определению автокорреляционной матрицы она Эрмитова и состоит из элементов АКФ (но я тут уже не помню, нужно почитать-помоделировать, забыл). SVD хорошо описаны у Голуба и Лоана в книжке (есть на русском кстати, могу поделиться)

спасибо, сейчас почитаю


именно о методе наименьших квадратов и идет речь.


буду весьма признательным, если поделитесь!

Если дело касается выравнивания и DFE, то там часто LDU разложение произведения тёплицевых матриц используется, поэтому A^H * A надо искать.

Вот книга (11 МБ djvu, не пролазит как прикрепленная)


Посмотрите у Хайкина главу про этот метод, если еще не читали. Там как раз описывается способ вычисления коэффициентов через псевдоинверсию матрицы данных. В отличие от некоторых "быстрых" методов инверсии обладает полной вычислительной устойчивостью и работает даже с сильно вырожденными матрицами (например в задачах с АФАР). Далее там есть глава про вычисления QR и SVD, но как раз у Голуба и Лоана более подробно это расписывается. Есть смысл покопать в гугле в сторону SVD Теплицовых матриц, наверняка есть какие то оптимизированные алгоритмы.


Ну тогда можно еще проще сделать: Ф = 1/N * sum_N(а*а^H), где а - выборка данных в момент времени t.

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

Если эти матрицы из "реальной жизни" , то есть шанс на то, что их собственные числа окажутся действительными. Оно так происходит с эрмитовыми. В таких случах может оказаться полезным метод Хаусхолдера, приводящий матрицу к тридиагональной форме. Он же, кстати, обычно является первой обязательной стадией QR/QL и SVD разложений, т.к. итерации делают не по плотной матрице, а по три- или дву- диагональным формам.

Есть задачи, в которых уже лишь приведение к 3- или 2-диагональным формам является вполне достаточным для получения эффективного (в смысле числа операций) решения, без необходимости доводить разложение до диагональной матрицы.

Xenia, для QR/QL/SVD декомпозиций можно использовать только вращение Гивенса, без Хаусхолдера. Или можно комбинировать Хаусхолдера и Гивенса для ускорения вычисления? Можете поделиться какими нибудь статьями с примерами?

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

Не прокатывает Там при поиске MSE решения для рассчета весов требуется разложить R = 1/snr*I + H*conj(H), H - теплицева (построена из ИХ канала). Т.е. в оценке входной корр. матрицы нет смысла - она в каком-то смысле известна (так как известна ИХ канала)

Хы, да это же похоже на MMSE для OFDM/MIMO OFDM , слагаемое H*conj(H) может быть найдено двояко, т.к. по формулировке MMSE это аппроксимация Rxx автокорреляционной матрицы, а она может быть найдена по алгоритму sample matrix inversion (никогда не понимал почему он так называется), как 1/N * sum_N(h*h^H). А если нет, то почему?

Не совсем похоже. В OFDM матрицы циркулярные (там халява), собственно за тем и придумали OFDM . Здесь они теплицевы - каждая H в формуле - это теплицева матрица.

PS всегда думал, что sample тут относится к тому, что корр матрица получена как оценка по входной выборке, типа из сэмплов.

Век живи, век учись). Я думал, что для OFDM Rxx = Rhh = E[h*h^H], где h - оценка канала в частотной области, полученная в пилотных точках. Вектор h мы берем на каждом OFDM символе. Далее w_MMSE = (Rxx + (1/SNR)*I)^-1 * Rxx, h_MMSE = h * w_MMSE^H. Ну вот как-то так. Мне не очень прозрачно как тут переписать матрицу Rhh через матрицу данных, т.к. в самом общем (и худшем) случае h произвольно дышит от блока к блоку.

А для MIMO получаем w_MMSE = (H*conj(H) + (1/SNR)*I)^-1 * conj(H), но тут H - матрица MIMO канала для данной поднесущей, которая в общем (и лучшем) случае произвольная комплексная матрица с линейно-независимыми строками.

ЗЫ: а про какой алгоритм вы говорите? Ниразу не попадалось инфы, чтобы MMSE применялся для SC систем.

Любую SC систему можно выравнивать эквалайзером в частотной области, если между блоками вставляется циклический префикс, нули или известная последовательность. В этом случае преобразование Фурье помогает в факторизации. Собственно поэтому это работает для OFDM. Я вобщем-то говорил про алгоритм оптимального формирования весов у DFE конечной длины. Познакомиться с ним можно у Al-Dhahir, например здесь и еще в куче статей

Насколько я знаю, метод Хаусхолдера применялся для этой цели еще со времен "Справочника" (Уилкинсон, Райнш). Оттуда же попал в большинство известных мне матпакетов, где для нахождения собственных значений и векторов приходилось вызывать не одну функцию, а последовательно три штуки: приведение к тридиагональной, QR-итерации и обратные итерации для восстановления соб.векторов. Но потом, когда в пакетах появилась для этих целей одна функция, стало сложно определить, как она внутри себя устроена. Например, MATLAB пишет, что qr() это его встроенная функция и кода не показывает. Тем не менее, подозреваю, что и там используется старый добрый Хаусхолдер. Вот и функция hess(), вычисляющая форму Гессенберга, использует все те же отражения Хаусхолдера, хотя его именем не называется.

По части применения в микроконтроллерах взгляните на это:
"Применение QR алгоритма для расчёта собственных структур корреляционных матриц на платформе ADSP-TS201 фирмы Analog Devices"
- там тоже Хаусхолдер.

И, наконец, цитата из "Алгебраическая проблема собственных значений – Часть 64":

Уилкинсон Дж.Х., Алгебраическая проблема собственных значений(1970):
что, видимо, объясняет причину популярности метода отражений.

Всегда смотрел в сторону хардверной платформы . Для реализации на FPGA/ASIC большую популярность сыскал именно метод Гивенса, т.к. он позволяет гибко параллелить и масштабировать вычисления в виде систолического массива элементарных вычислитилей (поворотных матриц 2х2). При этом в целочисленной арифметике можно сделать разложение без умножителей только на кордике.

В принципе вычисление по Гивенсу также переносятся и на процессорную архитектуру, но тогда элементарные повороты придется делать в цикле, вот тогда видимо и получается

С точки зрения арифметики операции в обоих методах подобны, значит и Хаусхолдера можно сделать без умножений. Интересно сохранится ли эта пропорция в таком случае?..

ЗЫ:

Когда разбирался со всей этой кухней написал простенькую реализацию qr по Гивенсу с плавучкой для верификации реализации с фиксированной точкой. Оказалось, что встроенная матлабовская qr даёт несколько большую ошибку для аргументов типа single чем самописная реализация. Матрицы комплексные и близкие к вырожденным в обоих случаях.

andyp, Xenia, спасибо за ссылки!

насколько я понял, лучше для оптимизации формулы X = (A^H * A)^-1 * A^H не отдельно оптимизировать умножение матриц и инверсию, а использовать псевдоинверсию.
но будет ли это менее затратно?

пусть матрица А имеет размерность [m x n], где m>n. используя псевдоинверсию и svd разложение получаем

X = V*S*U^H, где V - матрица порядка n, U - матрица порядка m, S - диагональная матрица размерностью [n x m];

перемножение этих матриц, плюс немалые затраты на svd разложение(с которым я еще до конца не разобрался)... будет ли это менее затратно, чем прямое вычисление по формуле X = (A^H * A)^-1 * A^H ?


кстати, матрица A у меня не совсем теплицевая.
теплицевыми являются матрицы A((1:2:end),: ) и A((2:2:end),: ).

но само перемножение A^H * A я могу делать пользуясь преимуществом теплицевой матрицы, т.к.
A^H * A = A((1:2:end),: )^H*A((1:2:end),: )+A((2:2:end),: )^H*A((2:2:end),: );

SVD операция затратная, но:
1) псевдоинверсия вычислительно устойчивый алгоритм
2) при переходе к матрице (A^H * A) вы будете работать с энергиями, т.к. это корреляционная матрица, разрядность вырастет в 2 раза, в общем случае вы не сможете откинуть добавочные разряды, т.к. потеряете информацию, которая в ней заключается (особенно это заметно в случае близком к вырождению, например АФАР, где собственные числа полезного сигнала много меньше собственных чисел помех).
3) SVD даст вам собственные числа матрицы данных, которые можно использовать с пользой в дальнейшем.
4) матрица данных у вас должна быть все таки Теплицевой или я чего то не понимаю, вы можете попробовать выиграть от SVD Теплицевой матрицы, а перейдя к (A^H * A) вы получете квадратную матрицу общего вида (она будет близка к Эрмитовой и тем ближе чем более длинные матрицы данных вы будете брать, но в общем случае она не Эрмитова, т.к. это только аппроксимация, поэтому я не уверен, что здесь можно выиграть от использования ее структуры)
5) если цель - алгоритм с фиксированной точкой, то по SVD есть литература как получить устойчивую имплементацию, хотя это будет сложно и весело. Можно в целых числах сделать инверсию (A^H * A) пользуясь QR декомпозицией, которая проще (и я бы даже мог поделиться матлабовскими скриптами по fixed point реализации), но см пункт 2 о росте разрядности.

В общем то, что я тут понаписал не претендует на истину в первой инстанции, т.к. я сам далеко не спец) но это мои соображения)

Странно у тебя получается. Вообще-то говоря, твоя X и есть псевдоинверсия матрицы A. Для ее вычисления нужно найти svd(A):

A = U*S*VH

S - диагональная

Тогда X = V*inv(S)*UH;

Инверсия S - просто инверсия элементов на главной диагонали.

andyp, мне кажется у него другая дилема: считать ли X = (A^H * A)^-1 * A^H напрямую или через X = V * S^-1 *U^H, с предварительным SVD. Мне кажется, что если бы вычисление напрямую было бы менее трудоёмким, то про 2й метод и не говорили бы

именно так


это радует

Она у меня не теплицевая, т.к. состовляется из двух разных импульсных характеристик. (на один символ приходится два отсчета.)
Если удалить из матрицы ряды через один, получиться теплицевая.


Интересует алгоритм с плавающей точкой

Вроде бы SVD будет наиболее эффективным вариантом, если требуется вычислить псевдоинверсию при этой структуре матрицы A. Кстати, такие "частично тёплицевы" матрицы вылазят в fractionally spaced эквалайзере.
Если это какой-то сорт MSE, то, может статься, что псевдоинверсия в замкнутом виде там и не нужна. Оптимальный вариант - обойтись без нее, что возвращает нас к статьям Al-Dhahir . У него там метод back substitution описан.

В моем случае используется DDE. Возможно ли в таком случае обойтись без псевдоинверсии? Может посоветуете какую-нибудь статью?

Так я ссылку приводил выше по теме именно на то, что нужно - быстрый способ вычисления весов в эквалайзере с обратной связью по решениям. Остальные роботы этого автора доступны здесь
http://www.utdallas.edu/~aldhahir/pubs.html

Виктор39, DDE = data directed equalizer? Который Nato STANAG*4285, Second example of demodulation technique?

да

UPD. Пардон, затупил. Произведение комплексных Теплицевых матриц вида A^H * A даёт Эрмитову, а вот при оценке автокорреляционной матрицы с увеличением времени оценки Эрмитова матрица стремится к реальной корреляционной у которой элементы на диагоналях, параллельных главной, равны (см. определение корреляционной матрицы).

В документе STANAG матрица весов W Теплицева, произведение W^H * W - Эрмитова матрица. Для обращения такой матрицы можно применить разложение Холецкого, которое легко сделать в арифметике с плавающей точкой. Разложение Холецкого устойчиво и требует меньше вычислений чем QR или SVD. Прямая инверсия матрицы тогда не понадобится, т.к. преобразование W^H * W к L*L^H приводит искомую систему линейных уравнений к двум треугольным. Про оптимизацию произведения Теплицевых матриц говорилось в начале темы.

объясните, почему инверсия матриц в таком случае не нужна?

Инверсия делается для решения системы линейных уравнений. Если систему можно переписать так, чтобы исключить прямую инверсию матрицы коэффициентов, то инверсия не нужна.
Удобным является приведение к треугольной форме, тогда система решается методом подстановки снизу/сверху. Разложение Холецкого представляет симметричную/Эрмитову матрицу в виде двух треугольных матриц: (W^H * W) = A = (L * L^H).
Исходная система: (W2^H * W2)*b = [W2^H * (r - W1*a)]
Полученная система: (L2 * L2^H)*b = [W2^H * (r - W1*a)]
Раскладываем её на 2 системы: L2 * t = [W2^H * (r - W1*a)] и L2^H * b = t.
И решаем их последовательно.

Общее замечание:
Воспоминания из студенческого прошлого говорят мне, что любую систему линейных уравнений можно атаковать например методом Гаусса. К инверсии можно обращаться, когда нужно решить дофига систем Ax = b для разных b. Матлаб, на сколько помню, тоже всегда советовал использовать x = A\b вместо x = inv(A)*b

Может статься, будет интересно
http://scicomp.stackexchange.com/questions...square-matrices

Познавательно насчёт матлаба


Это с точки зрения математики. А если смотреть на реализацию, то добавляются новые углы зрения. Для нас самый большой интерес представляют Теплицевы матрицы данных и Эрмитовы матрицы апроксимации автокорреляционной матрицы.

1) Если мы делаем обработку на процессоре с плавающей запятой, то можно решать СЛУ методом Холецкого, т.к. алгоритм имеет малое в сравнении с QR/LU для матриц общего вида количество вычислений; выполняется последовательно
2) Если мы делаем обработку на процессоре с фиксированной точкой, то возможно (?) выгоднее будет сделать QR/LU методом Гивенса или Хаусхолдера
3) Если мы делаем обработку на ПЛИС, то можно решать СЛУ через QR или SVD псевдоинверсию методом Гивенса, т.к. можно считать факторизацию параллельно на систолическом массиве из однообразных элементов
4) Еще есть класс быстрых алгоритмов, но они не все обладают устойчивостью и некоторые дают выигрыш только для очень больших матриц, поэтому это больше интересно для разгона научных вычислений

спасибо за помощь

не забываем про супербыстрый алгоритм Шура, описанный тут:
The Generalized Schur Algorithm for the Superfast Solution of Toeplitz Systems http://www.math.niu.edu/~ammar/papers/gsa.pdf