Hur kommer det sig att BLAS är optimerat?

Jag har hållit på med LAPACK den sista veckan. LAPACK är ett bibliotek för olika rutiner hur man kan beräkna massa saker på det allt möjliga konstiga sättet. Men det fungerar. LAPACK är ett uråldrigt bibliotek som uppdateras fortfarande. Det är skrivet i Fortran 90, men en utgåva har översats till ANSI C. Värt att notera är att det finns mycket C-kod där som är skrivet på K&R standarden

Hur som helst. I LAPACK så finns det blas t.ex. rutinen sgemm. Denna utför matrismultiplikation.
Jag öppnade sgemm.c filen och märkte att det var mycket kod för att beräkna C = A*B

Då är frågan. Hur kan denna fil vara optimerad om det är så mycket kod? Jag testade att klistra in koden i Godbolt och jag fick fram flera tusen rader assemblerkod. Denna kan väll inte vara optimerad? Jag menar, två stycken for-satser och man har gjort en matrismultiplikation. Enligt mig i alla fall.

Intressant! Jag gillar verkligen sådant. Jag kör LAPACK som ej tar hänsyn till hårdvara, dvs ingen GPU, Nvidia eller något annat sådant flashigt. Bara ren processorkraft. Jag har inte så mycket data att beräkna, typ 5000x5000 matriser.

Just nu har jag lite problem när jag kör ssyevd. Jag får nämligen felmeddelandet On entry to SSYEVD, parameter number 8 had an illegal value, vilket betyder att 8:e argumentet i funktionen ssyevd har ett felaktigt värde. (fixat)

Så åter till tråden.
Om jag kan köra en matrismultiplikation med två stycken for-satser. Är då SGEMM snabbare, om vi endast får använda vanlig processorkraft?

Ja, det är översatt. Jag skulle aldrig tro att någon skulle skriva en sådan oläsbar kod medvetet.

LAPACK har vissa rutiner i C som man inte hittar någon annanstans. Det finns boken "Numerical Recipies in C", men den innehåller bara det mest grundläggande. Vill man ha lite mera värre saker, så måste man använda LAPACK.

Man märker att koden är gammal om den innehåller K&R-kod.

Jag använder LAPACK.
Nu försöker jag hitta ett alternativ till att räkna ut SVD för osymmetriska matriser som är stora. Tidskrävande! Undra om jag ska använda rutinen för att beräkna egenvärden istället. Tydligen så är detta ett betydligt bättre sätt att räkna ut SVD än vanliga SVD.

function [U,S,V] = svdecon(X)
  % Input:
  % X : m x n matrix
  %
  % Output:
  % X = U*S*V'
  %
  % Description:
  % Does equivalent to svd(X,'econ') but faster
  %
  % Vipin Vijayan (2014)
  %X = bsxfun(@minus,X,mean(X,2));
  [m,n] = size(X);
  if  m <= n
      C = X*X';
      [U,D] = eig(C);
      clear C;

      [d,ix] = sort(abs(diag(D)),'descend');
      U = U(:,ix);

      if nargout > 2
          V = X'*U;
          s = sqrt(d);
          V = bsxfun(@(x,c)x./c, V, s');
          S = diag(s);
      end
  else
      C = X'*X;
      [V,D] = eig(C);
      clear C;

      [d,ix] = sort(abs(diag(D)),'descend');
      V = V(:,ix);

      U = X*V; % convert evecs from X'*X to X*X'. the evals are the same.
      %s = sqrt(sum(U.^2,1))';
      s = sqrt(d);
      U = bsxfun(@(x,c)x./c, U, s');
      S = diag(s);
  end
end

Nuvarande rutin är sgesvd och den är superseg. Har du ett bättre förslag, eller tycker du jag ska använda ssyevd istället?

Det tar 1163.852051 sekunder för mig att räkna ut en 5000x3000 matrix med sgesvd.

Antingen kan jag använda sgesdd som också räknar ut SVD fast på ett annat sätt. Vad tror ni är bäst för mig?
https://www.intel.com/content/www/us/en/docs/onemkl/code-samp...

Det är egenvektorerna som är målet. Det är PCA som är metoden. Orsaken varför PCA använder sig av SVD har med att PCA tar in en osymmetrisk matris, vilket eig i MATLAB ej kan räkna.

Men jag kommer inte köra direkt PCA, utan jag kommer köra Kernel PCA, vilket har visat sig vara en av dom bästa dimensionsreduceringsalgoritmerna då det finns valfrihet att utföra projektion samt använda olika kärnor beroende på datat.

I GNU octave så tog det 222 sekunder att göra special-SVD på en 5000x5000 matris. Jag vet inte vilken rutin EIG i GNU Octave använder.

X = randn(5000, 5000);
>> tic;[U, S, V] = svdecon(X);toc
Elapsed time is 222.039 seconds.

Med sgessd rutinen för en 5000x5000 matris tog det 440 sekunder. Det är självklart en förbättring jämfört med sgesvd som tog över 1000 sekunder för en 5000x3000 matris, inte ens fyrkantig!

I GNU Octave skulle jag aldrig kunna beräkna 5000x5000 matris med vanliga SVD.

Så undra om det blir bättre om man implementerar special SVD med LAPACK? Då kanske det går ännu snabbare?

Med ssyevd for en symmetrisk matris i C så kom jag till 217 sekunder för 5000x5000.

Jag tror inte jag ska testa special SVD i Lapack, dvs räkna ut egenvektorerna från en symmetrisk matris.
Detta har med att Lapack har stöd för att bara räkna ut U-matrisen från SVD, något som jag bara behöver.
Sedan så är multiplikation C = A*B väldigt krävande.

Jag testade jämföra sgemm, som utför C = alpha*A*B + beta*C, där alpha = 1 och beta = 0 för mig. Jämfört med två for-loopar.
sgemm visade sig kunna utföra multiplikation 150% snabbare än vanlig 2 for-satser.

Så det blir nog att hitta en svd-algoritm som passar mitt syfte.

Just nu har jag

• sgesdd - För generella matriser - Divide & Conquer algoritm - Har visat sig snabb

• sgesvd - För generella matriser - QR factorisering - Har visat sig seg, vet ej varför

• ssyevd - För symmetriska matriser - Divide & Conquer - Har visat sig snabb

Sedan finns så många olika algoritmer. Vem ska man ha, vet jag inte.

Jag implementerar BLAS rutiner nu. Men jag skulle behöva hjälp med att välja vilken SVD algoritm jag ska använda.

from time import time
import torch

# Quick and dirty implementation of Matlab/Octave tic/toc...
start = None

def tic():
    global start
    start = time()

def toc():
    global start
    elapsed = time() - start
    print(f'Elapsed time is {elapsed:.3f} seconds.')

tic()
A = torch.randn(5000, 5000)
U, S, V = torch.svd(A)
toc()

Elapsed time is 13.860 seconds.

Snyggt!
Men jag använder LAPACK. Eller menar du att MKL finns i C kod som ej tar hänsyn till hårdvara?

Enligt denna sida så finns det tre typer av SVD algoritmer för generella matriser. Jag tror inte det finns någon LAPACK algoritm för symmetriska SVD matriser då symmetriska SVD matriser och symmetriska egenvärdesmatriser är samma sak.

Men det som var kul att veta är att det finns tre typer av SVD algoritmer. Två av dessa är sgesdd, sgesvd. Tydligen så är, enligt Netlib, sgesdd bättre än sgesvd, oavsett läge.

https://netlib.org/lapack/lug/node71.html

Angående MKL....är den hårdvarubaserad?
Dvs MKL är förkompilerad, så man behöver bara länka. Men denna kod som är kompilerad, har väll tagit hänsyn till en specifik hårdvara?

Jag försöker hålla mig till ANSI C och koden ska kunna köras på alla plattformar. Så jag får nog svälja om det går segt för mig. Men så länge det går snabbare än mitt GNU Octave. Då har jag lyckats bra.

Det finns följande algoritmer:
sgejsv, sgesdd, sgesvd, sgesvdq, sgesvdx

Dessa kan utföra SVD för generella matriser. Ryktet på nätet säger att sgesdd är alltid bättre än sgesvd. Detta håller jag med om också. Rutinen sgejsv är bara till för precision. Då finns det två algoritmer kvar, nämligen sgesvdq, sgesvdx.

Det som skiljer sgesvdq, sgesvdx från sgesdd är att sgesdd räknar både ut U, S, V, eller S endast. Medan sgesvdq och sgesvdx så har man flera valfriheter att välja vad man vill räkna ut. Jag är bara intresserad utav U matrisen i mitt PCA-fall.

https://netlib.org/lapack/explore-html/d4/dca/group__real_g_e...

Vilken algoritm bör jag använda: sgesvdq, sgesvdx ? Eller ska jag hålla mig till sgesdd?

Edit:
Glöm det! I mitt Lapack 3.2.1 så finns inte sgesvdq, sgesvdx tillgängliga.

Fungerar MKL på alla hårdvaror, oavsett ålder?

Men om man inte kör x86 då? Man kanske kör AVR eller något annat mystiskt? Jag håller på med ett projekt som handlar om att porta numerisk linjär algebra till plattformsoberoende system.

Jag märker att den LAPACK kod jag använder, är mycket väl skriven på 80-talet. Vissa delar i alla fall.

Ja, det var detta jag ville ha svar på. BLAS och LAPACK är alltså optimerad för specifik hårdvara. I detta fall använder jag BLAS och LAPACK som är plattformsoberoende. Det är därför det går lite segt. Men med lite tur så kan jag komma snabbare än mitt GNU Octave, vilket är det primära målet.

Jag vet dock inte om mitt GNU Octave är så väl optimerat. Det gick ju segt igår att räkna ut SVD.

@Elgot

Nu har jag hittat lite mer information.
För beräkning av PCA så fungerar sgesvd bättre än sgesdd, trots att sgesdd är snabbare om man ska ha både U, S, V. Men vill man bara ha U, så är sgesvd snabbare då man har valet att välja om man vill bara beräkna U.

Här är ett litet program som jag har gjort som beräknar PCA

>> tic; [P, W] = mi.pca(X, 2);toc
The size of the matrix is 5000x5000. Do you want to apply PCA onto it? 1 = Yes, 0 = No: 1
Elapsed time is 237.496 seconds.
>>

X är en 5000x5000 matris och jag vill dimensionera ned X till 2 dimensioner. Ca 240 sekunder. Rätt dåligt. Men jag tror jag använder vanliga BLAS i mitt Octave.

För sgesvd med 5000x5000 matris så tog det 256 sekunder.
För sgesdd med 5000x5000 matris så tog det 349 sekunder.

Men jag tror jag är dålig på att välja optimering.