WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

  
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Value Function Iteration mit GPU Parallelisierung lösen
 

Kaja

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 28.11.2019, 15:09     Titel: Value Function Iteration mit GPU Parallelisierung lösen
  Antworten mit Zitat      
Hallo zusammen

Ich würde gerne eine simple Matlab Funktion schreiben, welche Value Function Iteration
parallelisiert auf dem GPU berechnen kann.
Ich möchte, ungefähr, diesen Code für Julia auf Matlab replizieren:
https://github.com/giob1994/GPU-Ope.....under%20Uncertainty.ipynb

Dabei habe ich besonders Problem zu verstehen, wie man diesen Teil des Codes auf Matlab implementiert:

Code:
# Write kernel for GPU manually:
gpu_call(grid, (grid, V, policy, z, P, Float32(alpha), Float32(beta), Float32(delta), Float32(sigma), UInt32(SIZE_GRID), UInt32(SIZE_Z)))
do state, grid, V, policy, z, P, alpha, beta, delta, sigma, SIZE_GRID, SIZE_Z
# Each kernel executes for one value of the capital grid:
idx = @linearidx grid


Gibt es ähnliche Funktionen wie "gpu_call" und "@linearidx" für Matlab?
Die parallel computing toolbox habe ich installiert : )

Die ähnlichste Funktion, welche ich finden konnte, ist
Code:
KERN = parallel.gpu.CUDAKernel(PTXFILE,CPROTO)

jedoch bezieht sich diese Funktion so weit ich verstehe auf C for CUDA oder OpenCL CODE, welcher ich nicht verstehe.

Dies ist ungefähr, wie die Funktion schlussendlich aussehen sollte:
Code:
function [V,pol] = VFI_own_gpu_attempt(alpha,beta,delta,eta,z_grid,k_grid,pi_z,tol)
size_k_grid = size(k_grid,1);
size_z_grid = length(z_grid);

k_grid_G = gpuArray(k_grid);
z_grid_G = gpuArray(z_grid);
pi_z_G = gpuArray(pi_z);
V0 = ones(size_k_grid,size_z_grid,'gpuArrays');
V = ones(size_k_grid,size_z_grid,'gpuArrays');
pol = zeros(size_k_grid,size_z_grid,'gpuArrays');

while abs(V-V0)>tol
V0 = V;
% write kernel
%gpu_call(...)

% each kernel executes for one value of the capital grid
%idx = @linearidx grid

for i_z = 1:size_z_grid
F = -Inf;
pol_i = uint(1)
    for i_k = 1_size_k_grid
    c = z_grid_G(i_z)*k_grid_G(idx)^alpha + (1-delta)*k_grid_G(idx) - k_grid_G(i_k)M
        if c>0
        F0 = ((c)^(1-eta)-1)/(1-eta)
            for j = 1:size_z_grid
                F1 = F0 + beta*pi_z_G(i_z,j)*V(i_k,j);
            end
        end
        if F1 > F
        F = F1;
        pol_i = uint64(i_k);
        end
    end
V(idx,i_z) = F;
pol(idx,i_z) = pol_i;

end


Ich bedanke mich im Voraus für jegliche Hilfe! Smile


Harald
Forum-Meister

Forum-Meister

Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 28.11.2019, 15:23     Titel:
  Antworten mit Zitat      
Hallo,

für elementweise Operationen auf der GPU würde ich arrayfun empfehlen.
Dazu muss man nur die Daten anfangs auf die GPU bringen, wie du das ja schon gemacht hast. Der Rest geht "automatisch" auf der GPU.

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Kaja

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 28.11.2019, 16:02     Titel:
  Antworten mit Zitat      
Hallo Harald

Danke vielmals für die rasche Antwort.

Ich bin mir nicht ganz sicher, wie ich arrayfun hier richtig brauchen würde.

Wäre dies so richtig:
1. Im main.m vor der VFI-Funktion:
Code:

k_grid_G = gpuArray(k_grid);
z_grid_G = gpuArray(z_grid);
...

2. Die Funktion rufen als arrayfun:
[V,policy]=arrayfun(VFI_own_attempt, V0, V, alpha,beta,delta,eta,z_grid,k_grid,pi_z,tol)
3. Die VFI_own_attempt.m Funktion wäre dann die selbe wie ohne GPU, alles Zuweisungen zum GPU ('gpuArray') geschehen ausserhalb der Funktion.

Was ich nicht genau verstehe:
Bei der Dokumentation für arrayfun werden Beispiele für die Funktion immer mit syms gezeigt.

Code:
B = arrayfun(@(x), func(x,_,_)  ,__)


Welche Funktion würden diese syms für meine VFI-Funktion einnehmen?

Besten Dank nochmals! Smile
 
Harald
Forum-Meister

Forum-Meister

Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 28.11.2019, 16:16     Titel:
  Antworten mit Zitat      
Hallo,

genau, alles in der äußeren for-Schleife muss in die Funktion für arrayfun.

arrayfun hat an sich nichts mit syms zu tun, und ich finde auch keine Seite, die das kombiniert. Die Seite, nach der du dich richten solltest, ist diese hier:
Code:


Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Kaja

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 28.11.2019, 17:38     Titel:
  Antworten mit Zitat      
Danke vielmals! Very Happy

Ich habe das jetzt einmal versucht so zu implementieren.
Leider stoppt Matlab wenn die Linie mit arrayfun aufgerufen wird mit:
Zitat:
Error using gpuArray/arrayfun
Matrix dimensions must agree.

Error in own_gpu_attempt (line 36)
[V,pol]= arrayfun(@VFI_own_gpu_attempt1, V0, alpha, beta, delta, eta,size_k, size_z,z,k,pi_z_G)


Ich sehe nicht wirklich, wass ich hier falsch tue:

Code:
alpha=0.35;
beta= 0.984;
delta=0.01;
eta=2;

tol=10^(-4);

% Not using Tauchen but some transition matrix with 3 grid points for stochastic process
z_grid = [0.4,0.8,1.2];
pi_z = [0.7,0.2,0.1;0.1,0.8,0.1;0.05,0.1,0.85];

% capital grid
k_min = 0.01;
k_max = 20;
n_k = 100;
k_grid = linspace(k_min,k_max,n_k);

% prepare objects for GPU
k = gpuArray(k_grid);
z = gpuArray(z_grid);
pi_z_G = gpuArray(pi_z);

size_k = size(k,1);
size_z = length(z);

V0 = zeros(size_k,size_z,'gpuArray');
V = ones(size_k,size_z,'gpuArray');
pol = zeros(size_k,size_z,'gpuArray');
diff=max(max(V-V0));

while diff>tol
V0 = V;
[V,pol]= arrayfun(@VFI_own_gpu_attempt1, V0, alpha, beta, delta, eta,size_k, size_z,z,k,pi_z_G)
diff=max(max(V-V0));
end


und die Funktion

Code:
function [V,pol] = VFI_own_gpu_attempt1(V0,alpha, beta, delta, eta,size_k,size_z,z,k,pi_z_G)

y_k = (k.^alpha)'*z+ (1-delta)*k(ones(1,size_z),:)';

for i_k = 1:size_k
    for i_z = 1:size_z

        low_k=1;
     
    if(y_k(i_k,i_z) > k(end))
        high_k = length(k);
    else
        high_k= find(k > y_k(i_k,i_z), 1);
    end
   
    if(k(high_k) > y_k(i_k,i_z))
   high_k = high_k -1;
    end

    N_k = high_k;

%maximization
F = ((y_k(i_k,i_z)*ones(N_k,1) - k(low_k:high_k)').^(1-eta))/(1-eta) + beta*V0(low_k:high_k,:)*pi_z_G(i_z,:)';
[V(i_k,i_z),pol(i_k,i_z)]=max(F);


end

end


Hast du eine Idee, wo der Fehler herkommen könnte?
Siehst du noch andere Probleme in diesem Code?

Liebe Grüsse
 
Harald
Forum-Meister

Forum-Meister

Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 28.11.2019, 17:43     Titel:
  Antworten mit Zitat      
Hallo,

die Schleifen und auch das y_k müssen auf jeden Fall aus der Funktion raus. Das komponentenweise Anwenden übernimmt ja arrayfun. y_k sollte Input sein, auf dem elementweise gearbeitet wird.

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Kaja

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 28.11.2019, 18:25     Titel:
  Antworten mit Zitat      
Hm, danke für den Hinweis.

Ich bekommen immernoch den selben Fehlerhinweis.

Den Code habe ich so angepasst:
Code:
while diff>tol
V0 = V;
y_k = (k.^alpha)'*z+ (1-delta)*k(ones(1,size_z),:)';

for i_k = 1:size_k
    for i_z = 1:size_z
        low_k=1;
     
    if(y_k(i_k,i_z) > k(end))
        high_k = length(k);
    else
        high_k= find(k > y_k(i_k,i_z), 1);
    end
   
    if(k(high_k) > y_k(i_k,i_z))
   high_k = high_k -1;
    end

    N_k = high_k;

[V,pol]= arrayfun(@VFI_own_gpu_attempt1, V0,y_k,i_k,i_z,high_k,low_k, alpha, beta, delta, eta,size_k, size_z,z,k,pi_z_G)
diff=max(max(V-V0));
    end
end
end

und
Code:
function [V,pol] = VFI_own_gpu_attempt1(V0,y_k,i_k,i_z,high_k,low_k,alpha, beta, delta, eta,size_k,size_z,z,k,pi_z_G)

%maximization
F = ((y_k(i_k,i_z)*ones(N_k,1) - k(low_k:high_k)').^(1-eta))/(1-eta) + beta*V0(low_k:high_k,:)*pi_z_G(i_z,:)';
[V(i_k,i_z),pol(i_k,i_z)]=max(F);


end
 
Harald
Forum-Meister

Forum-Meister

Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 28.11.2019, 19:22     Titel:
  Antworten mit Zitat      
Hallo,

sorry, das war wohl unklar: die Schleifen über i_k und i_z müssen komplett weg. Der if-Teil muss dagegen in die Funktion rein.

Meine Empfehlung wäre, den Versuch mit arrayfun zunächst auf der CPU zu testen. Dort hast du bessere Debugging-Möglichkeiten.

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Kaja

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 28.11.2019, 22:35     Titel:
  Antworten mit Zitat      
Danke für die Tipps!

Auf dem CPU habe ichs zum Laufen gebracht (ist aber irgendwo noch ein Fehler).
Wenn ich aber diesen CPU Code:
Code:

diff=max(max(V-V0));

while abs(diff)>tol
V0 = V;
y_k = (k.^alpha)'*z+ (1-delta)*k(ones(1,size_z),:)';

for i_k = 1:size_k
    for i_z = 1:size_z
             
[V,pol]= VFI_own_gpu_attempt1(V0,y_k,i_k,i_z, alpha, beta, delta, eta,size_k, size_z,z,k,pi_z_G);
    end
end
diff=max(max(V-V0));
end
 

Zum GPU Code wechsle, kommt wieder dieselbe Fehlermeldung.

Code:
while abs(diff)>tol
V0 = V;
y_k = (k.^alpha)'*z+ (1-delta)*k(ones(1,size_z),:)';

for i_k = 1:size_k
    for i_z = 1:size_z
             
[V,pol]= arrayfun(@VFI_own_gpu_attempt1,V0,y_k,i_k,i_z, alpha, beta, delta, eta,size_k, size_z,z,k,pi_z_G);
    end
end
diff=max(max(V-V0));
end


Die Funktion sieht immernoch so aus:

Code:
function [V,pol] = VFI_own_gpu_attempt1(V0,y_k,i_k,i_z,alpha, beta, delta, eta,size_k,size_z,z,k,pi_z_G)
    low_k=1;
    if(y_k(i_k,i_z) > k(end))
        high_k = length(k);
    else
        high_k= find(k > y_k(i_k,i_z), 1);
    end
   
    if(k(high_k) > y_k(i_k,i_z))
   high_k = high_k -1;
    end

    N_k = high_k;
%maximization
F = ((y_k(i_k,i_z)*ones(N_k,1) - k(low_k:high_k)').^(1-eta))/(1-eta) + beta*V0(low_k:high_k,:)*pi_z_G(i_z,:)';
[V(i_k,i_z),pol(i_k,i_z)]=max(F);


end


Ich verstehe nicht wirklich wie die Matrixdimensionen im CPU-Fall aufgehen können aber im GPU-Fall nicht. Scheint mathematisch nicht wirklich aufzugehen.

Liebe Grüsse
 
Harald
Forum-Meister

Forum-Meister

Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 28.11.2019, 23:08     Titel:
  Antworten mit Zitat      
Hallo,

der Punkt ist, arrayfun auf der CPU zu verwenden und das zu debuggen.

Und bitte die Ratschläge beachten:
Zitat:
die Schleifen über i_k und i_z müssen komplett weg

(wenn du arrayfun verwenden willst)

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Kaja

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 28.11.2019, 23:13     Titel:
  Antworten mit Zitat      
Hallo Harald

Entschuldigung, da habe ich dich wirklich falsch verstanden. (Ich dachte, du meinst "weg aus der Funktion".)

Meinst du mit "arrayfun auf CPU ausführen":
- die Funktion gleich behalten
- die Objekte (k_grid,z_grid, etc.) **nicht** als normales double lassen?

Wie kann ich die beiden for-schleifen entfernen?
Mit Vektorization?

Liebe Grüsse
 
Kaja

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 28.11.2019, 23:21     Titel:
  Antworten mit Zitat      
das **nicht** sollte natürlich weg... sorry
 
Harald
Forum-Meister

Forum-Meister

Beiträge: 24.448
Anmeldedatum: 26.03.09
Wohnort: Nähe München
Version: ab 2017b
     Beitrag Verfasst am: 28.11.2019, 23:27     Titel:
  Antworten mit Zitat      
Hallo,

Zitat:
Wie kann ich die beiden for-schleifen entfernen?

Das ist doch genau der Punkt an arrayfun, dass du dann eben keine Schleifen mehr brauchst, weil arrayfun die Funktion automatisch elementweise ausführt.

Grüße,
Harald
_________________

1.) Ask MATLAB Documentation
2.) Search gomatlab.de, google.de or MATLAB Answers
3.) Ask Technical Support of MathWorks
4.) Go mad, your problem is unsolvable ;)
Private Nachricht senden Benutzer-Profile anzeigen
 
Kaja

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 29.11.2019, 00:55     Titel:
  Antworten mit Zitat      
Ja, habe realisiert, dass ich nicht verstehe, wie arrayfun funktioniert. Dachte es sei eine auf GPU spezialisiert Funktion.

Konnte mit Hilfe der Dokumentation nicht wirklich verstehen, wie **arrayfun** korrekt zu gebrauchen ist. Es gibt sehr wenig kurze Beispielcodes oder Erklärvideos zu dieser Funktion.

Ich informier mich morgen noch einmal mehr.
 
Kaja

Gast

Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
     Beitrag Verfasst am: 29.11.2019, 13:48     Titel:
  Antworten mit Zitat      
Hallo Harald

Ich habe den Code jetzt so gut wie möglich abgeändert:

Code:
tic
alpha=0.35;
beta= 0.984;
delta=0.01;
eta=2;

tol=10^(-4);

% Not using Tauchen but some transition matrix with 3 grid points for stochastic process
z_grid = [0.4,0.8,1.2];
pi_z = [0.7,0.2,0.1;0.1,0.8,0.1;0.05,0.1,0.85];

% capital grid
k_min = 0.01;
k_max = 20;
n_k = 1000;
k_grid = linspace(k_min,k_max,n_k);

% prepare objects for GPU
% k = gpuArray(k_grid);
% z = gpuArray(z_grid);
% pi_z_G = gpuArray(pi_z);
k = k_grid;
z = z_grid;
pi_z_G = pi_z;

size_k = size(k,2);
size_z = length(z);

V0 = zeros(size_k,size_z);
V = ones(size_k,size_z);
pol = zeros(size_k,size_z);
diff=max(max(V-V0));

while abs(diff)>tol
V0 = V;
y_k = (k.^alpha)'*z+ (1-delta)*k(ones(1,size_z),:)';

             
[V,pol]= arrayfun(@VFI_own_gpu_attempt1,V0,y_k,alpha, beta, delta, eta,size_k, size_z,z,k,pi_z_G);
diff=max(max(V-V0));
end
toc


Und die Funktion:

Code:
function [V,pol] = VFI_own_gpu_attempt1(V0,y_k,alpha, beta, delta, eta,size_k,size_z,z,k,pi_z_G)
    low_k=1;
    if(y_k > k(end))
        high_k = length(k);
    else
        high_k= find(k > y_k, 1);
    end
   
    if(k(high_k) > y_k)
   high_k = high_k -1;
    end

    N_k = high_k;
%maximization
F = ((y_k*ones(N_k,1) - k(low_k:high_k)').^(1-eta))/(1-eta) + beta*V0(low_k:high_k,:)*pi_z_G';
[V,pol]=max(F);


end


Der Fehlerhinweis lautet:

Code:
Error using arrayfun
All of the input arguments must be of the same size and shape.
Previous inputs had size 1000 in dimension 1. Input #4 has size 1

Error in own_gpu_attempt (line 42)
[V,pol]= arrayfun(@VFI_own_gpu_attempt1,V0,y_k,alpha, beta, delta, eta,size_k, size_z,z,k,pi_z_G);


Das Problem scheint, dass bei arrayfun alle inputs dieselbe Dimension haben müssen.

Also schreibe ich die Funktion so um:
Code:
function [V,pol] = VFI_own_gpu_attempt1(V0,y_k,z,k,pi_z_G)
    size_k = size(k,2);
    size_z = length(z);
    alpha=0.35;
    beta= 0.984;
    delta=0.01;
    eta=2;
   
    low_k=1;
    if(y_k > k(end))
        high_k = length(k);
    else
        high_k= find(k > y_k, 1);
    end
   
    if(k(high_k) > y_k)
   high_k = high_k -1;
    end

    N_k = high_k;
%maximization
F = ((y_k*ones(N_k,1) - k(low_k:high_k)').^(1-eta))/(1-eta) + beta*V0(low_k:high_k,:)*pi_z_G';
[V,pol]=max(F);


end


Gibt es einen Trick, wie ich nun pi_z_G, z und k auch in die arrayfun einspeisen kann?
Irgendwie mit z=repmat(z,1000)?
Bin mir nicht sicher ob dies funktioniert und mathematisch aufgeht.

Liebe Grüsse
 
Neues Thema eröffnen Neue Antwort erstellen

Gehe zu Seite 1, 2  Weiter
Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

 Du kannst Beiträge in dieses Forum schreiben.
Du kannst auf Beiträge in diesem Forum antworten.
Du kannst deine Beiträge in diesem Forum nicht bearbeiten.
Du kannst deine Beiträge in diesem Forum nicht löschen.
Du kannst an Umfragen in diesem Forum nicht mitmachen.
Du kannst Dateien in diesem Forum posten
Du kannst Dateien in diesem Forum herunterladen
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.