Версия для печати темы Нажмите сюда для просмотра этой темы в оригинальном формате

Форум программистов > Алгоритмы > blas и lapack в python

Автор: mrgloom 23.6.2014, 09:12

Допустим имеем обёртку вокруг BLAS на python (т.е. имеем numpy), можем ли руками написать ф-ии из LAPACK?(по минимуму мне надо eig, svd.) я знаю про numpy.linalg.svd и numpy.linalg.eig так же там можно посмотреть сорцы, которые показывают в случае svd Цитата     The decomposition is performed using LAPACK routine _gesdd , а в случае eig Цитата     implemented using the _geev LAPACK routines но там правда есть ссылка на Цитата     References     ----------     G. Strang, *Linear Algebra and Its Applications*, 2nd Ed., Orlando, FL,     Academic Press, Inc., 1980, Various pp. Для чего всё это нужно: переписать стандартные функции из numpy (svd и eig) с использованием numpy.memmap и cudamat (тот же blas обертка вокруг cublas как я понимаю). Так же для поддержки больших матриц можно использовать блочное перемножение матриц (хотя может еще какие то ф-ии blas придется переделать? и я чего то не учёл?) Что точнее нужно: может быть есть lapack собранный только из вызовов blas? (так ли это и возможно ли это?) возможно есть примеры на других языках? (т.к. нет никакого желания разбираться в фортран коде и ф-иях с названиями типа dgemm) какие книги по линейной алгебре(для инженеров скорее, а не для математиков) лучше читать, которые более близки к реальности и рассказывают о подводных камнях(например вырожденные случаи (типа сингулярная матрица) или там потеря точности (типа таким методом никто не считает, а вот более точный) и т.д.), опять же учитывая, что blas и обёртки у нас уже есть и требуется написать только высокоуровневые алгоритмы.

Автор: _Y_ 24.6.2014, 19:27

Даже интересно: а при чем здесь ветка "Алгоритмы"?

Автор: mrgloom 25.6.2014, 08:42

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

Автор: Pavia 25.6.2014, 21:54

mrgloom,  Если вопрос практический, то лучше его было задать в разделе питона.  По поводу фортрана. Библиотеке LAPACk и BLAS есть написанные на Си. Если ручками, то  честно даже не знаю что вам подсказать.  Цитата может быть есть lapack собранный только из вызовов blas? (так ли это и возможно ли это?)  А из чего он ещё состоит?!!! Цитата(mrgloom @  25.6.2014,  08:42 ) а во-вторых нужна рекомендация книги по практической линейной алгебре, которая как бы тоже про алгоритмы. Из книг трудно подсказать могу вот автора подсказать Дж. Голуб. Если найдёшь свисни. Цитата(mrgloom @  25.6.2014,  08:42 ) ну во-первых тут алгоритмы svd и eig и может кто то поделиться информацией, если их реализовывал руками. Проще стенку головой пробить. Чем разобраться как их сделать.  По началу кажется просто. А потом выясняется что  описание идёт с использованием терминов высшей математики. Приходиться её учит. Затем что процесс только итерационный QR. Как следствие бороться с устойчивостью. Что для симметричный матриц надо симметричность поддерживать. Что для действительных матриц оказываются решения бывают комплексные. А вот с ними я так и не разобрался. Толи есть там алгоритм, то ли его нет.  А ещё у матриц при умножении их есть проблема с переполнением чисел Float даже на маленькой матрице 8х8  Поэтому все и пользуются готовым . А в готовых используется другой алгоритм нежели чем тот который описан в книгах.  Так что помимо QR ещё и процесс(аглоритм) Грамма-Шмитта надо смотреть.

Автор: Фантом 25.6.2014, 23:52

Цитата(mrgloom @  23.6.2014,  10:12 ) Допустим имеем обёртку вокруг BLAS на python (т.е. имеем numpy), можем ли руками написать ф-ии из LAPACK?(по минимуму мне надо eig, svd.)  Что значит "можем"? Python в принципе позволяет это сделать, сможете ли Вы - это другой вопрос... Цитата(mrgloom @  23.6.2014,  10:12 ) может быть есть lapack собранный только из вызовов blas? (так ли это и возможно ли это?) Нет, и это невозможно. Цитата(Pavia @  25.6.2014,  22:54 ) По поводу фортрана. Библиотеке LAPACk и BLAS есть написанные на Си. Это тоже обертки исходного фортрановского кода. Цитата(Pavia @  25.6.2014,  22:54 ) Проще стенку головой пробить. Чем разобраться как их сделать.  Это не настолько уж страшно, но все-таки здравое зерно в этом утверждении есть. При необходимости в написании софта такого рода куда проще (и полезнее для дела) разобраться в Фортране, чем пытаться наперекор судьбе чесать правой ногой левое ухо и реализовывать вычислительные алгоритмы на языке, совершенно не предназначенном для подобных задач. Я уж не говорю о том, что даже если почесаться удастся, то работающий результат из-за всех промежуточных оберток будет иметь производительность раз так в десять меньшую, чем нормальный код.

Автор: Pavia 26.6.2014, 05:28

Цитата(Фантом @  25.6.2014,  23:52 ) Это тоже обертки исходного фортрановского кода Нет не обёртка, а с транслированный код.  Цитата(Фантом @  25.6.2014,  23:52 )  разобраться в Фортране, чем пытаться наперекор судьбе чесать правой ногой левое ухо и реализовывать вычислительные алгоритмы на языке, совершенно не предназначенном для подобных задач.  Это вы просто программировать не умете. Язык этот как раз хорошо подходит для подобных задач.  Цитата(Фантом @  25.6.2014,  23:52 ) Я уж не говорю о том, что даже если почесаться удастся, то работающий результат из-за всех промежуточных оберток будет иметь производительность раз так в десять меньшую, чем нормальный код. Нет обёртки влияют на производительность как 0.01-0.001 раз и то обычно даже ещё меньше. Нет никаких 10 раз.  А вот Python не быстрый язык и писать на нём я бы не стал.

Автор: mrgloom 26.6.2014, 09:13

Цитата(Pavia @  25.6.2014,  21:54 ) Цитата может быть есть lapack собранный только из вызовов blas? (так ли это и возможно ли это?)  А из чего он ещё состоит?!!! Цитата(Фантом @  25.6.2014,  23:52 ) Цитата(mrgloom @  23.6.2014,  10:12 Найти цитируемый пост) может быть есть lapack собранный только из вызовов blas? (так ли это и возможно ли это?) Нет, и это невозможно. ну так вот я и не понял, по сути blas это набор разных частных случаев перемножения матриц? что еще присутствует в lapack? какие то дополнительные ф-ии хэлперы? т.е. я не понимаю достаточно ли ф-ии из blas для использования их как минимальных неделимых кирпичиков для построение стандартных алгоритмов линейной алгебры?  код действительно есть транслированный на С, но в нем так же сложно разобраться(ф-ии с 15 аргументами и т.д.) можно использовать прямой вызов  http://docs.scipy.org/doc/scipy-dev/reference/linalg.blas.html http://wiki.scipy.org/PerformanceTips единственный способ похоже это поковыряться в самих ф-иях lapack https://groups.google.com/forum/#!topic/julia-dev/mmgO65i6-fA кстати пишут, что  Цитата The routine dgesdd uses a divide and conquer algorithm for the bidiagonal SVD, whereas dgesvd uses a QR algorithm. In summary, I recommend always calling the dgesdd if the memory is available, otherwise dgesvd.

Автор: Фантом 26.6.2014, 11:25

Цитата(Pavia @  26.6.2014,  06:28 ) Нет не обёртка, а с транслированный код.  Вы исходники CBLAS когда-нибудь видели? Если нет - я покажу. Вот, например, cblas_dtbmv.c (случайно ткнул в первый попавшийся файл), начальные строки: Код /*  * cblas_dtbmv.c  * The program is a C interface to dtbmv.  *  * Keita Teranishi  5/20/98  *  */ #include "cblas.h" #include "cblas_f77.h" void cblas_dtbmv(const enum CBLAS_ORDER order, const enum CBLAS_UPLO Uplo,                  const enum CBLAS_TRANSPOSE TransA, const enum CBLAS_DIAG Diag,                  const int N, const int K, const double  *A, const int lda,                  double  *X, const int incX) {    char TA;    char UL;    char DI; #ifdef F77_CHAR    F77_CHAR F77_TA, F77_UL, F77_DI; #else    #define F77_TA &TA    #define F77_UL &UL    #define F77_DI &DI    #endif #ifdef F77_INT    F77_INT F77_N=N, F77_lda=lda, F77_K=K, F77_incX=incX; #else    #define F77_N N    #define F77_K K    #define F77_lda lda    #define F77_incX incX #endif    extern int CBLAS_CallFromC;    extern int RowMajorStrg;    RowMajorStrg = 0;    CBLAS_CallFromC = 1;    if (order == CblasColMajor)    {       if (Uplo == CblasUpper) UL = 'U';       else if (Uplo == CblasLower) UL = 'L';       else        {          cblas_xerbla(2, "cblas_dtbmv","Illegal Uplo setting, %d\n", Uplo);          CBLAS_CallFromC = 0;          RowMajorStrg = 0;          return;       }       if (TransA == CblasNoTrans) TA = 'N'; Как, похоже на самостоятельный код? Или все-таки на обертку? Цитата(Pavia @  26.6.2014,  06:28 ) Это вы просто программировать не умете. Язык этот как раз хорошо подходит для подобных задач.     Цитата(Pavia @  26.6.2014,  06:28 ) А вот Python не быстрый язык и писать на нём я бы не стал.  Юноша, Вы перед декларированием громких заявлений хотя бы для себя определитесь, подходит Python для вычислительных задач или нет. А то сами себе противоречите.

Автор: Фантом 26.6.2014, 12:04

Цитата(mrgloom @  26.6.2014,  10:13 ) ну так вот я и не понял, по сути blas это набор разных частных случаев перемножения матриц? Практически да (еще сложения).  Цитата(mrgloom @  26.6.2014,  10:13 ) что еще присутствует в lapack? какие то дополнительные ф-ии хэлперы?  Решение задач линейной алгебры - это и решение линейных систем (тут одними умножениями и сложениями не обойтись), решение задачи о собственных числах матрицы, задачи факторизации и т.д. и т.п. Цитата(mrgloom @  26.6.2014,  10:13 ) код действительно есть транслированный на С, но в нем так же сложно разобраться(ф-ии с 15 аргументами и т.д.) Поскольку в природе есть f2c, то формально получить код C несложно, но он будет, во-первых, почти нечитаем, во-вторых, неэффективен. Честно же писать все это заново на C желающих нет - долго, сложно и никому не нужно.

Автор: mrgloom 25.7.2014, 09:49

Появился вопрос как сравнить работу 2-х SVD алгоритмов, т.е. один точный стандартный, но медленный ,его берем за основу ,второй допустим использует некие аппроксимации и поэтому менее точный, но более быстрый.  1.Можно сравнить собственные векторы и значения: Тут показана разность полученных результатов: Цитата u               [,1]          [,2]          [,3]          [,4]         [,5]  [1,] 2.775558e-16  2.220446e-16  1.387779e-15 -7.771561e-16  0.037587204  [2,] 0.000000e+00  1.387779e-16 -1.165734e-15  2.026157e-15  0.062777989  [3,] 1.665335e-16  1.387779e-17 -1.949829e-15 -1.165734e-15 -0.034826470  [4,] 0.000000e+00 -8.326673e-17  4.440892e-16  2.130240e-15 -0.027517156  [5,] 2.220446e-16  4.440892e-16  5.551115e-16  1.998401e-15  0.006668889  [6,] 1.110223e-16 -5.551115e-17 -6.938894e-17  6.106227e-16 -0.008904566  [7,] 5.551115e-17 -1.110223e-16  5.551115e-16  1.477984e-15 -0.006729478  [8,] 1.110223e-16 -1.249001e-16 -5.273559e-16 -2.359224e-16 -0.143823784  [9,] 1.665335e-16 -1.665335e-16 -5.828671e-16  7.771561e-16 -0.043560560 [10,] 1.665335e-16  0.000000e+00 -1.498801e-15 -1.665335e-16 -0.008980835  v               [,1]          [,2]          [,3]          [,4]        [,5] [1,]  3.330669e-16 -2.220446e-16 -1.942890e-16  5.551115e-16  0.11045048 [2,]  8.881784e-16  1.332268e-15 -5.551115e-17 -2.706169e-16  0.09073229 [3,]  0.000000e+00  0.000000e+00 -1.942890e-16  1.221245e-15 -0.05183135 [4,]  5.551115e-17  6.661338e-16  1.498801e-15 -1.831868e-15  0.03964778 [5,]  1.665335e-16  4.996004e-16 -3.219647e-15 -3.108624e-15 -0.06788847 [6,]  2.220446e-16  1.942890e-16  2.359224e-15  4.274359e-15  0.01698281 [7,]  4.440892e-16  1.054712e-15 -4.440892e-16 -1.387779e-16  0.10417172 [8,] -1.110223e-16 -1.026956e-15  3.330669e-16 -3.330669e-16 -0.09256007   d [1]  0.000000e+00  1.110223e-15  1.332268e-15 -2.220446e-16  6.660773e-03 2. Можно сравнить ошибку реконструкции изначальной матрицыhttp://en.wikipedia.org/wiki/Matrix_norm#Frobenius_norm 3. Можно проверить вектора на ортогональность. Но вопрос в том, какая точность\разность в выдаваемых значениях приемлема для реальных задач?

Powered by Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)