Более эффективный способ вычисления XSX ^ T и XSy

Question

Мне нужно вычислить следующие матрицы: M = XSX ^ T а также V = XSy

я хотел бы знать, что более эффективная реализация использует blas, зная, что S является симметричной и определенной положительной матрицей измерения n, X имеет m строк и n столбцов, а y - вектор длины n.

Моя реализация следующая: Я вычисляю A = XS с использованием dsymm, а затем с помощью dgemm получаем M = AX ^ T, в то время как dgemv используется для получения V = Ay.

Я думаю, что, по крайней мере, M можно вычислить более эффективным способом, поскольку я знаю, что M симметрично и определенно положительно.

Answer 1

Не зная деталей вашей проблемы, может быть полезно распознать нули в матрицах. Разделение матриц для достижения этого может обеспечить значительные преимущества. Является ли M суммой многих подматриц XSX?

Для V = XSy, где y - вектор, а X и S - матрицы, вычисление S.y тогда X. (Sy) должно быть лучше, если X.S не является необходимым вычислением для M.

Answer 2

Ваш код - лучшее, что BLAS может сделать для вас. Нет операции BLAS, которая может использовать тот факт, что M симметрична.

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

Могу ли я узнать о размерах? И позвольте мне также вдохновить некоторые другие источники для повышения производительности: Собственная сборка вашей реализации BLAS, сравнение с MKL, ACML, OpenBLAS, ATLAS. Очевидно, вы можете написать свою собственную версию, которая будет использовать AVX, FMA intrinsics. Вы должны быть в состоянии сделать лучше, чем некоторые обобщенные библиотеки. Кроме того, какова точность вашей переменной с плавающей точкой?

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

Но вы не ответили на мой вопрос относительно реализации BLAS и типа машины.

Опять же, оптимизация за пределами этого пункта невозможна без дополнительных навыков :(. А если серьезно, не беспокойтесь об этом. Есть причина, почему BLAS не делает оптимизацию, которую вы запрашиваете. Стоит хлопот. Пойдите со своим решением.

И не забывайте инвестировать в использование поплавков, а не удваивать. На R конвертируй все, чтобы плавать. Для команд Lapack используйте только sgemX