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

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

{ 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)

Сообщение отредактировал andyp - Apr 10 2015, 09:28

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

Сообщение отредактировал serjj - Apr 10 2015, 09:48

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


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


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

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

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


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


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

Сообщение отредактировал serjj - Apr 10 2015, 11:35

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

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

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

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

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

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

Сообщение отредактировал andyp - Apr 10 2015, 13:52

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

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

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

Сообщение отредактировал andyp - Apr 10 2015, 14:07

Век живи, век учись). Я думал, что для 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, спасибо за ссылки!

Сообщение отредактировал serjj - Apr 11 2015, 08:08