Я ищу наиболее эффективный код IDL для замены оператора умножения матрицы IDL (#) для конкретной диагонально-ориентированной (не диагональной или диагонально-симметричной) матрицы с 3 различными значениями: единица на диагонали; единица плюс дельта справа от диагонали; единица минус такая же дельта влево.
Проблемная область
IDL (фиксированный; не подлежит обсуждению; извините); смазывание изображения на беззатворной ПЗС-системе формирования изображения.
Основная постановка проблемы
Дано
матрица 1024x1024 "EMatrix" с единицей по диагонали; (1-дельта) слева от диагонали; (1+дельта) вправо; дельта = 0,044.
другая матрица 1024x1024, Изображение
Запрос
- какой самый быстрый код IDL для расчета (изображение # EMatrix)?
2014-09-16: см. обновление ниже
Задний план
Более крупная постановка задачи (из которой умножение матрицы кажется только самой медленной частью, и оптимизация всей процедуры не повредит):
Раздел 9.3.1.2 (страница PDF 47; внутренняя страница 34) и другие документы в том же каталоге (извините, как новичок, я могу опубликовать только две ссылки)
Моя работа до сих пор
Это сейчас (2014-09-26) примерно на порядок быстрее, чем оператор IDL # для матрицы 1024x1024.
Подробности
Наивная операция — O(n^3) и выполняет около миллиарда (2^30) умножений с двойной точностью и примерно такое же количество сложений; Далее Википедия сообщает мне, что алгоритм Штрассена сокращает это значение до O (n ^ 2,807), или ~ 282M+ умножений для n = 1024.
Разбивая это для простого случая 3x3, скажем, изображение и EMatrix
image EMatrix
[ 0 1 2 ] [ 1 p p ]
[ 3 4 5 ] # [ m 1 p ]
[ 6 7 8 ] [ m m 1 ]
где p представляет собой 1+дельта (1,044), а m представляет собой 1-дельта (0,956).
Из-за повторения m и p должно быть доступно упрощение: глядя на средний столбец для изображения, результат для трех строк должен быть
[1,4,7] . [1,p,p] = m*(0) + 1*1 + p*(4+7)
[1,4,7] . [m,1,p] = m*(1) + 1*4 + p*(7)
[1,4,7] . [m,m,1] = m*(1+4) + 1*7 + p*(0)
куда . представляет точечный (внутренний?) продукт, т.е. [1,4,7].[a,b,c] = (1a + 4b + 7c)
Исходя из этого, вот что я сделал до сих пор:
Средний член — это просто сам столбец, а суммы, умноженные на m и p, очень похожи на кумулятивные суммы смежных участков столбца, возможно, перевернутых (для m), сдвинутых на единицу и с первым элементом, установленным на ноль.
то есть для м:
; shift image down one:
imgminusShift01 = shift(image,0,1)
; zero the top row:
imgminusShift01[*,0] = 0
; Make a cumulative sum:
imgminusCum = total( imageShift01, 2, /cumulative)
Для p imgplusShift01Cum следует по существу тому же пути, но с поворотом (..., 7) до и после, чтобы переворачивать элементы вверх и вниз.
Когда у нас есть эти три матрицы (исходное изображение с его дочерними элементами imgminusShift01Cum и imgplusShift01Cum), желаемый результат будет
m * imgminusShift01Cum + 1 * image + p * imgplusShift01Cum
Для выполнения сдвигов и поворотов я использую индексы, которые сдвигаются и поворачиваются сами по себе и сохраняются в общем блоке для последующих вызовов, что экономит еще 10-20%.
Резюме, пока
2014-09-16: см. обновление ниже
И это дает ускорение 5+.
Я ожидал немного большего, потому что я думаю, что умножил 3M и добавил 2M, так что, возможно, выделение памяти — дорогостоящая часть, и я должен делать что-то вроде REPLICATE_INPLACE (мой IDL ржавый — я мало что делал с 6.4 и ранний 7.0).
Или, может быть, матричное умножение IDL лучше теории?
Альтернативные подходы и другие мысли:
Можем ли мы что-то сделать с тем фактом, что EMatrix равна единице плюс матрица с нулями по диагонали и +/- дельта в верхнем и нижнем треугольниках?
Суммируя по столбцам, я последовательно обращаюсь к данным «неправильным» способом, но действительно ли это сэкономит время, если сначала перенести изображение (это всего 8 МБ)?
Очевидно, что выбор другого языка, получение массива графических процессоров или написание DLM были бы другими вариантами, но давайте пока оставим это в IDL.
advСПАСИБО (да, я такой старый ;-),
Брайан Карчич 20 июля 2014 г.
Обновление 2014-09-16
Подход Диего привел меня к гораздо более простому решению:
Я думаю, мы поняли; мой первый проход был слишком сложным, со всеми вращениями.
Используя обозначения Диего, но транспонированные, я ищу
[K] = [IMG] # [E]
Поскольку это умножает столбцы [IMG] на строки [E], между столбцами [IMG] нет взаимодействия, поэтому для анализа нам нужно посмотреть только один столбец [IMG], скалярные (внутренние) произведения которого , со строками [E], становятся одним столбцом результата [K]. Расширяя эту идею до одного столбца матрицы 3x3 с элементами x, y и z:
[Kx] [x] [1 1+d 1+d]
[Ky] = [y] # [1-d 1 1+d]
[Kz] = [z] [1-d 1-d 1 ]
Рассмотрим конкретно элемент Ky выше (один элемент [K], соответствующий y в [IMG], разбитый на формулу, использующую только скаляры):
Ky = x * (1-d) + y * 1 + z * (1+d)
Обобщение для любого y в столбце любой длины:
Ky = sumx * (1-d) + y * 1 + sumz * (1+d)
Где скаляры sumx и sumz представляют собой суммы всех значений выше и ниже y соответственно в этом столбце [IMG]. Н.Б. sumx и sumz — это значения, относящиеся к элементу y.
Перестановка:
Ky = (sumx + y + sumz) - d * (sumx - sumz)
Ky = tot - d * (sumx - sumz)
куда
tot = (sumx + y + sumz)
т. е. tot — это сумма всех значений в столбце (например, в IDL: tot = total(IMG,2)).
Так что до сих пор я в основном дублировал работу Диего; остальная часть этого анализа преобразует это последнее уравнение для Ky в форму, пригодную для быстрой оценки в IDL.
Решение уравнения tot для sumz:
sumz = tot - (y + sumx)
Замена обратно в Ky:
Ky = tot - (sumx - (tot - (y + sumx)))
Ky = tot - ((2 * sumx) + y - tot)
Ky = tot + (tot - ((2 * sumx) + y)
Использование sumxy для представления суммы всех значений в столбце сверху вниз до y включительно (IDL: [SUMXY] = total([IMG],2,/CUMULATIVE))
sumxy = sumx + y
а также
sumx = sumxy - y
Замена обратно в Ky:
Ky = tot + (tot - ((2 * (sumxy - y)) + y)
Ky = tot + (tot + y - (2 * sumxy))
Таким образом, если мы можем оценить tot и sumxy для каждого элемента [IMG], т.е. если мы можем оценить матрицы [TOT] и [SUMXY], и у нас уже есть [IMG] как матричная версия y, тогда это простая линейная комбинация этих матриц.
В IDL это просто:
[SUMXY] = TOTAL([IMG],2,/CUMULATIVE)
[TOT] = [SUMXY][*,N-1] # REPLICATE(1D0,1,N)
т.е. [TOT] — это последняя строка [SUMXY], продублированная для формирования матрицы из N строк.
И окончательный код выглядит так:
function lorri_ematrix_multiply,NxN,EMatrix
NROWS = (size(NxN,/DIM))[1]
SUMXY = TOTAL(NxN,2,/CUMULATIVE)
TOT = SUMXY[*,NROWS-1] # REPLICATE(1,NROWS,1d0)
RETURN, TOT + ((EMatrix[1,0] - 1d0) * (TOT + NxN - (2d0 * SUMXY)))
end
который в нашей системе просто на порядок быстрее, чем [IMG]#[E].
Н.Б. дельта = (EMatrix[1,0] - 1d0)
Ву ху!
