[c++ vs matlab] ottimizzare codice c++ per configgere matlab

[mamo139]

Ciao a tutti!
spero che l'argomento possa essere di interesse a qualcuno visto che per la mia testolina da programmatore della domenica non è banale :zizi:

allora sto scrivendo un software che attraverso simulazione di montecarlo trovi il prezzo di un certo derivato esotico (precisamente una opzione magrabe ma non credo vi sia di interesse)

l'ho fatto sia in matlab che in c++, sicuro che il c++ avrebbe stracciato matlab di tanto... e invece mi ritrovo con un codice matlab piu veloce

eccovi le due funzioni incriminate

c++

//using normal distribution object from Technical Report 1
typedef std::tr1::ranlux64_base_01 Myeng; 
typedef std::tr1::normal_distribution<double> ndistr;

double Margrabe_option_MC_price(double r, double t, double s1, double s2, double q1, double q2, 
								double std1, double std2, double corr, double iterations){
	//variables and objects initialization
	register long x;
	double price = 0; //the final price will be stored here
	double payoff;
	double z1, z2; //to store generated random numbers
	double matr21, matr22;
	double sigma1, sigma2;
	double part1, part2, sqrtt;

	clock_t start, stop;

	Myeng eng;
	ndistr stdnorm(0, 1); //creating the standard normal distribution generator

	//seeding the random number generator
	eng.seed((unsigned long) time(NULL)); 
	
	//calculating sigmas
	matr21 = corr * std2;
	matr22 = sqrt(1-corr*corr) * std2;

	//precalculating fixed part of drive formula
	part1 = exp((r - q1 - 0.5*std1*std1)*t);
	part2 = exp((r - q2 - 0.5*std2*std2)*t);
	sqrtt = sqrt(t);

	for(x=0; x < iterations ;x++){
		//generating random numbers
		z1 = stdnorm(eng);
		z2 = stdnorm(eng);

		//creating sigmas
		sigma1 = z1 * std1;
		sigma2 = matr21 * z1 + matr22 * z2;

		//calculating payoff
		payoff = (s1 * part1 * exp(sqrtt * sigma1)) - (s2 * part2 * exp(sqrtt * sigma2));
		price += (payoff > 0) ?  payoff : 0;
	}
	price /= iterations;
	price *= exp((-r)*t);

	return price;
}

matlab

function [p,ci] = ExchangeMC(V0,U0,sigmaV,sigmaU,rho,T,r,NRepl)
eps1 = randn(1,NRepl);
eps2 = rho*eps1 + sqrt(1-rho^2)*randn(1,NRepl);
VT = V0*exp((r - 0.5*sigmaV^2)*T + sigmaV*sqrt(T)*eps1);
UT = U0*exp((r - 0.5*sigmaU^2)*T + sigmaU*sqrt(T)*eps2);
DiscPayoff = exp(-r*T)*max(VT-UT, 0);
[p,s,ci] = normfit(DiscPayoff);

ebbene matlab per fare 100.000 generazioni ci mette 24 centesimi di secondo, contro i 38 del codice c++.

Considerando che ho sempre pensato il c/c++ essere i linguaggi di alto livello piu performanti che esistono (correggetemi se sbaglio) devo aver sbagliato io qualcosa.

come rendereste voi il codice c piu efficiente?

da ignorante mi è venuto da pensare che il problema potrebbe essere nella funziona di generazione di numeri casuali di una normale standard!!
ho controllato e matlab lo risolve con questi calcoli
eps1 = randn(1,NRepl);
eps2 = rho*eps1 + sqrt(1-rho^2)*randn(1,NRepl);
in un elegante 8 centesimi

mentre la funzione che ho trovato io nelle librerie del boost che da quel che ho capito riscreano il pacchetto aggiuntivo Technical Report 1
for(x=0; x < iterations ;x++){
//generating random numbers
z1 = stdnorm(eng);
z2 = stdnorm(eng);
}
se ne esce con un goffo 26 centesimi di secondo

come posso migliorare questo aspetto?
inoltre avete altre osservazioni sui codici??

ci tengo particolarmente a rendere il c++ piu performante...

grazie
Marco

[shinya]

Originariamente inviato da mamo139

//using normal distribution object from Technical Report 1
typedef std::tr1::ranlux64_base_01 Myeng;

Non sono un gran esperto di numeri casuali, ma hai già provato anche con altri engine?

[starfred]

ciao, come suppongo tu abbia un pc multicore, potresti cercare di adottare una soluzione di calcolo parallelo multithreading.
Per esempio nel tuo codice il calcolo di z1 e sigma1 può essere messo in parallelo con il calcolo di z2.
Più parallelizzi e più hai possibilità che il programma vada più veloce.
(Ho detto possibilità in quanto 10 thread su un dual core non vanno quasi sicuramente meglio di 2 thread )

[mamo139]

Originariamente inviato da shinya

Non sono un gran esperto di numeri casuali, ma hai già provato anche con altri engine?

Ho provato un po di altri engine a caso ma nulla batte sostanzialmente questo qui...
qualcuno conosce qualche motore particolarmente efficiente?

ho il sospetto che matlab sia efficiente nel generare n numeri casuali contemporaneamente, piuttosto che n volte singolarmente

Originariamente inviato da starfred

ciao, come suppongo tu abbia un pc multicore, potresti cercare di adottare una soluzione di calcolo parallelo multithreading.
Per esempio nel tuo codice il calcolo di z1 e sigma1 può essere messo in parallelo con il calcolo di z2.
Più parallelizzi e più hai possibilità che il programma vada più veloce.
(Ho detto possibilità in quanto 10 thread su un dual core non vanno quasi sicuramente meglio di 2 thread )

Sinceramente ci avevo già pensato, mettere in parallelo x di queste funzioni dove x è il numero di core della macchina, tutte quante calcolando n/x iterazioni invece che n e poi unendo i risultati.

ma questa sarà l'ultima chicca, voglio almeno pareggiare prima di andare al multithreading

[shinya]

Originariamente inviato da mamo139

Ho provato un po di altri engine a caso ma nulla batte sostanzialmente questo qui...
qualcuno conosce qualche motore particolarmente efficiente?

Mersenne Twister?

[starfred]

se matlab usa il multithreading sarà molto difficile andar più veloce con un single thread... non credo che il C++ basti, potresti provare soluzioni asm miste c++

[mamo139]

Originariamente inviato da shinya

Mersenne Twister?

DUE MIGLIORIE:
ieri notte ho separato la generazione dei numeri random immagazzinandoli in un array, quindi ora ci sono due for al posto di uno ma c'è stata una miglioria:

da 0.44 a 0.34

questa mattina ho sostituito il tr1 con:
una versione del Mersenne Twister per generare numeri random e una funzione che pare essere molto ottimizzata per normalizzarli (spero sia fatta abbastanza bene da non crearmi problemi a livello di qualità dei dati )
per esequire 10 simulazioni da 1.000.000 ci mette ora 0.30 al posto di 0.34 del tr1
ma bisogna migliorare ancora

eccovi il codice allo stato attuale:

inline double Margrabe_option_MC_price(double r, double t, double s1, double s2, double q1, double q2, 
								double std1, double std2, double corr, long iterations){
	//variables and objects initialization
	register long x,y;
	double price = 0; //the final price will be stored here
	double payoff;
	double *z1; //to store generated random numbers
	double matr21, matr22;
	double sigma1, sigma2;
	double part1, part2, sqrtt;

	//memory allocation
	z1 = new double[iterations*2];

	//seeding the random number generator
#ifdef RANDOM_MODE_TR1
	Myeng eng;
	ndistr stdnorm(0, 1); //creating the standard normal distribution generator
	eng.seed((unsigned long) time(NULL));
#endif
#ifdef RANDOM_MODE_Mersenne_twister	
    init_genrand((int)time(NULL));
#endif
	//calculating sigmas
	matr21 = corr * std2;
	matr22 = sqrt(1-corr*corr) * std2;

	//precalculating fixed part of drive formula
	part1 = exp((r - q1 - 0.5*std1*std1)*t);
	part2 = exp((r - q2 - 0.5*std2*std2)*t);
	sqrtt = sqrt(t);

	for(x=0; x < iterations ;x++){
		y=x+x;

		//generating random numbers
		#ifdef RANDOM_MODE_TR1
		z1[y] = stdnorm(eng);
		z1[y+1] = stdnorm(eng);
		#endif
		#ifdef RANDOM_MODE_Mersenne_twister	
		z1[y] = normsinv(genrand_real3());
		z1[y+1] = normsinv(genrand_real3());
		#endif
	}

	for(x=0; x < iterations ;x++){
		y=x+x;

		//creating sigmas
		sigma1 = z1[y] * std1;
		sigma2 = matr21 * z1[y] + matr22 * z1[y+1];

		//calculating payoff
		payoff = (s1 * part1 * exp(sqrtt * sigma1)) - (s2 * part2 * exp(sqrtt * sigma2));
		if(payoff > 0) price += payoff;
	}
	price /= iterations;
	price *= exp((-r)*t);

	delete [] z1;

	return price;
}

Originariamente inviato da starfred

se matlab usa il multithreading sarà molto difficile andar più veloce con un single thread... non credo che il C++ basti, potresti provare soluzioni asm miste c++

si matlab di default usa multithreading in tutte le funzioni matematiche con vettori superiori a qualche migliaio di dati. Ma lo sto facendo andare su un 2 core, quindi possiamo batterlo secondo me

io purtroppo asm non lo conosco, ma se tu hai voglia di farmi vedere come potrei sostituire qualche pezzetto con codice asm apprezzerei molto (e poi così mi metto anche a vedere come funziona l'assembly)

[starfred]

ciao, dunque per prima cosa devi definire una funzione, per esempio creo una funzione che faccia y=x+x;
la chiamo somma che dovrà fare queste semplici cose
int somma(int x){
int y;
y=x+x;
return y;
}

Per prima cosa, per tradurra in asm bisogna fare la traduzione dei nomi; poiché tu usi C++ per comodità (anche se si potrebbe fare a meno) bisogna dichiararla extern quindi avremo
extern "C" somma(int x)

A questo punto passiamo all'asm e quindi al file .s che noi chiameremo fun.s
quindi

.text
.global _somma
_somma:

pushl %ebp
movl %esp, %ebp
subl $4, %esp

movl 8(%ebp), %eax
addl 8(%ebp), %eax


leave
ret

Ora, niente panico e facciamo chiarezza. Per prima cosa di semplice intuizione il .global serve per il nome e la successiva istruzione definisce la funzione.

Poi, salvo il contenuto del registro ebp, esso va salvato perché non vogliamo che il suo contenuto si perda. Successivamente sposto l'indirizzo base del record di attivazione che è in esp in ebp (dato che non vogliamo perdere l'indirizzo della pila che mi serve per "trovare" le variabili passate alla funzione).
A questo punto creo spazio per le variabili locali ovvero un int cioè 32 bit cioè 4 byte e quindi sposto il puntatore della pila di una posizione in "su" appunto per creare lo spazio per la variabile y.
A questo punto il mio record di attivazione è

d (valore) -4
<- punto si trova l'indirizzo EBP
EIP (+4)
x (valore che è stato passato alla funzione) (+8)

Quindi per prendere il valore x devo "scendere" nel mio record a +8 (ricordati sempre che c'è sempre EIP che occupa 4byte!!!) quindi ecco spiegato la movl 8(%ebp), %eax
ovvero "sposta il CONTENUTO della locazione di memoria identificato dall'INDIRIZZO contenuto in ebp incrementato di 8 byte in EAX".
Successivamente faccio banalmente la somma e con leave elimino lo spazio delle variabili locali ripristinando il valore di esp e successivamente di ebp.
Infatti il codice della leave equivale a

movl %ebp, %esp
popl %ebp

[-MiStO-]

probabilmente nel tuo caso è trascurabile ma io eviterei nei cicli for il post incremento(i++) a favore del preincremento(++i)
questo perchè il post crea un temp con il valore dell'indice originale, incrementa l'indice e ritorna il temp: hai quindi una variabile in più inutilizzata per ogni iterazione
inoltre in alcuni casi può essere utile procedere all'unroll dei cicli for, ma non so se questo è il caso (si può sempre tentare )

[marco.r]

ancora prima direi che vale la pena di dare una profilata per vedere dove viene speso il tmepo, senza tentare ottimizzazioni a casaccio

[WarDuck]

Probabilmente Matlab avrà alcune parti scritte in assembly e ottimizzate per l'uso di alcune estensioni dei processori moderni per il calcolo vettoriale.

[das]

Se scopri che la maggior parte del tempo viene spesa nella generazione dei numeri casuali prova ad usare una pennina come questa questo:


http://www.entropykey.co.uk/

Hai anche il vantaggio di avere numeri realmente casuali e non pseudocasuali come quelli ottenibili con metodi software.
Nel sito dice che incrementa moltissimo le prestazioni rispetto agli algoritmi software normalmente utilizzati.

[clockover]

Originariamente inviato da das

Se scopri che la maggior parte del tempo viene spesa nella generazione dei numeri casuali prova ad usare una pennina come questa questo:


http://www.entropykey.co.uk/

Hai anche il vantaggio di avere numeri realmente casuali e non pseudocasuali come quelli ottenibili con metodi software.
Nel sito dice che incrementa moltissimo le prestazioni rispetto agli algoritmi software normalmente utilizzati.

Si ma il tempo di I/O è di sicuro maggiore... bella discussione quella di strizzare per bene il codice per una gara con matlab

[-MiStO-]

Originariamente inviato da marco.r

ancora prima direi che vale la pena di dare una profilata per vedere dove viene speso il tmepo, senza tentare ottimizzazioni a casaccio

...e in effetti hai ragione
quindi mamo139 profila(vs2010, non ricordo se professional o ultimate, ha un discreto tool di profilazione integrato) e poi si vedrà dove mettere mano

[clockover]

Scusami ma perchè hai questa cosa

for(x=0; x < iterations ;x++){
		y=x+x;

		//generating random numbers
		#ifdef RANDOM_MODE_TR1
		z1[y] = stdnorm(eng);
		z1[y+1] = stdnorm(eng);
		#endif
		#ifdef RANDOM_MODE_Mersenne_twister	
		z1[y] = normsinv(genrand_real3());
		z1[y+1] = normsinv(genrand_real3());
		#endif
	}

	for(x=0; x < iterations ;x++){
		y=x+x;

		//creating sigmas
		sigma1 = z1[y] * std1;
		sigma2 = matr21 * z1[y] + matr22 * z1[y+1];

		//calculating payoff
		payoff = (s1 * part1 * exp(sqrtt * sigma1)) - (s2 * part2 * exp(sqrtt * sigma2));
		if(payoff > 0) price += payoff;
	}

non riesci a mettere tutto in un ciclo? A occhio sembra di si comunque... poi potrebbe sfuggirmi qualcosa

[mamo139]

alla fine ho risolto il problema e battuto matlab ottimizzato per multithreading (utilizzato su un dual core) con un codice C non ancora ottimizzato
ho trovato un super algoritmo di generazione e l'ho implementato! in questo caso matlab ha un ottimo algoritmo di generazione di numeri casuali!!

ecco l'algoritmo utilizzato per chiunque fosse interessato
http://www.mathworks.it/company/news...g01_cleve.html

avrò provato una decina di algoritmi trovati su dei paper prima di arrivare a questo e nessuno andava così bene! quindi se vi interessa l'argomento consiglio la lettura!

ora ottimizzo anche il mio codice C per andare in multithreading e vi posto i risultati e il codice finale!

------------------------------------------------------------------------------------

figa la chiavetta generatrice ma non ho fondi ne tempo per provarla

Originariamente inviato da clockover

Scusami ma perchè hai questa cosa

for(x=0; x < iterations ;x++){
		y=x+x;

		//generating random numbers
		#ifdef RANDOM_MODE_TR1
		z1[y] = stdnorm(eng);
		z1[y+1] = stdnorm(eng);
		#endif
		#ifdef RANDOM_MODE_Mersenne_twister	
		z1[y] = normsinv(genrand_real3());
		z1[y+1] = normsinv(genrand_real3());
		#endif
	}

	for(x=0; x < iterations ;x++){
		y=x+x;

		//creating sigmas
		sigma1 = z1[y] * std1;
		sigma2 = matr21 * z1[y] + matr22 * z1[y+1];

		//calculating payoff
		payoff = (s1 * part1 * exp(sqrtt * sigma1)) - (s2 * part2 * exp(sqrtt * sigma2));
		if(payoff > 0) price += payoff;
	}

non riesci a mettere tutto in un ciclo? A occhio sembra di si comunque... poi potrebbe sfuggirmi qualcosa

si si puo unire in un unico ciclo ma ho separato per poter profilare e vedere quale parte di calcoli consumava piu tempo
una profilazione un po casereccia , comunque la pardita di velocità che ho registrato da questa scorporazione è davvero trascurabile

Originariamente inviato da Antonio23

dove l'hai letta sta cosa?

sul sito ufficiale di matlab. l'ho anche verificata avviando matlab in modalità senza multithreading e ho verificato che le funzioni di operazioni fra vettori traggono grandissimi benefici dal multithreading
invece ad esempio rand no perchè non è ottimizzata (c'era scritto anche questo sul sito ufficiale, ma ho verificato ed è vero)

Originariamente inviato da -MiStO-

inoltre in alcuni casi può essere utile procedere all'unroll dei cicli for, ma non so se questo è il caso (si può sempre tentare )

scusami l'ignoranza, cosa intendi per unroll dei cicli for?

[Torav]

Originariamente inviato da mamo139

alla fine ho risolto il problema e battuto matlab ottimizzato per multithreading (utilizzato su un dual core) con un codice C non ancora ottimizzato
ho trovato un super algoritmo di generazione e l'ho implementato! in questo caso matlab ha un ottimo algoritmo di generazione di numeri casuali!!

ecco l'algoritmo utilizzato per chiunque fosse interessato
http://www.mathworks.it/company/news...g01_cleve.html

avrò provato una decina di algoritmi trovati su dei paper prima di arrivare a questo e nessuno andava così bene! quindi se vi interessa l'argomento consiglio la lettura!

ora ottimizzo anche il mio codice C per andare in multithreading e vi posto i risultati e il codice finale!

------------------------------------------------------------------------------------

figa la chiavetta generatrice ma non ho fondi ne tempo per provarla


si si puo unire in un unico ciclo ma ho separato per poter profilare e vedere quale parte di calcoli consumava piu tempo
una profilazione un po casereccia , comunque la pardita di velocità che ho registrato da questa scorporazione è davvero trascurabile



sul sito ufficiale di matlab. l'ho anche verificata avviando matlab in modalità senza multithreading e ho verificato che le funzioni di operazioni fra vettori traggono grandissimi benefici dal multithreading
invece ad esempio rand no perchè non è ottimizzata (c'era scritto anche questo sul sito ufficiale, ma ho verificato ed è vero)



scusami l'ignoranza, cosa intendi per unroll dei cicli for?

Io ti consiglierei di profilare il codice per sapere dove vale la pena ottimizzare. Io uso valgrind (con cachegrind) che ti dice quanto tempo spendi su ogni riga di codice. Veramente una bomba.
Probabilmente perderai parecchio tempo a fare esponenziali. Compili con -ffast-math? Se gli argomenti degli exp hanno un limite minore e maggiore puoi ottimizzare con una look up table. Ma prima ti conviene capire dov'è che il tuo programma passa la maggior parte del tempo.

[mamo139]

Originariamente inviato da Antonio23

quando intendi avviare matlab "in modalità mulththreading" vuol dire che avvii i WORKERS (in numero pari in genere a quanti core disponi?)? in pratica sfrutti le funzionalità del parallel toolbox? se è così quello non è matlab puro... il parallel toolbox si paga a parte...

neanche sapevo dell'esistenza di questo toolbox

comunque matlab implementa il multithreading in automatico per la maggior parte delle funzioni base su matrici. Se si vuole disabilitare il multithreading bisogna farlo all'avvio come opzione.

ma tanto lo batto con un core solo

MATLAB versions 7.8 and later require that you decide between threading and no threading at the time you launch your MATLAB session. Multithreading is enabled by default. To disable this feature, start MATLAB using the new singleCompThread option.

fonte: http://www.mathworks.it/help/techdoc/rn/bry1ecg-1.html

As of MATLAB 7.4 (R2007a), MATLAB supports multithreaded computation for a number of linear algebra functions (e.g. matrix multiply), element-wise numerical functions (e.g. cos), and expressions that are combinations of element-wise functions (e.g. y=4*x*(sin(x) + x^3)). These functions automatically execute on multiple threads and you do not need to explicitly specify commands to create threads in your code.



Element Wise Functions and Expressions:
------------------------------------------------------------------------------------------------
Functions that speed up for double arrays > 20k elements

1) Trigonometric: ACOS(x), ACOSH(x), ASIN(x), ASINH(x), ATAN(x), ATAND(x), ATANH(x), COS(x), COSH(x), SIN(x), SINH(x), TAN(x), TANH(x)

2) Exponential: EXP(x), POW2(x), SQRT(x)

3) Operators: x.^y
For Example: 3*x.^3+2*x.^2+4*x +6, sqrt(tan(x).*sin(x).*3+8);

Functions that speed up for double arrays > 200k elements

4) Trigonometric: HYPOT(x,y), TAND(x)

5) Complex: ABS(x)

6) Rounding and remainder: UNWRAP(x), CEIL(x), FIX(x), FLOOR(x), MOD(x,N), ROUND(x)

7) Basic and array operations: LOGICAL(X), ISINF(X), ISNAN(X), INT8(X), INT16(X), INT32(X)

Linear Algebra Functions:
------------------------------------------------------------------------------------------------
Functions that speed up for double arrays > 40k elements (200 square)


1)Operators: X*Y (Matrix Multiply), X^N (Matrix Power)

2)Reduction Operations : MAX and MIN (Three Input), PROD, SUM

3) Matrix Analysis: DET(X), RCOND(X), HESS(X), EXPM(X)

4) Linear Equations: INV(X), LSCOV(X,x), LINSOLVE(X,Y), A\b (backslash)

5) Matrix Factorizations: LU(X), QR(X) for sparse matrix inputs

6) Other Operations: FFT and IFFT of multiple columns of data, FFTN, IFFTN, SORT, BSXFUN, GAMMA, GAMMALN, ERF,ERFC,ERFCX,ERFINV,ERFCINV, FILTER

fonte: http://www.mathworks.it/support/solu...ution=1-4PG4AN