Numerics Libraries (Foundations)

Cuchulainn · November 29th, 2011, 7:11 pm

QuoteOriginally posted by: PolterCuch, "CholeskyNov2011.zip" file uses solely the LU decomposition. There is also (unreferenced, so didn't look at it) procedure "Cholesky".Any chance you could post the version that uses the Cholesky decomposition A = L * L^T = U^T * U ?Yes, thanks. I noticed that. Accidentally removed it when cleaning.Will add Cholesky and go through the same steps again. Now UD --> transpose(L). What I plan to do here is to code the 10 most common numerical linear algebra routines in uBLAS and compare with other libraries and integrate into applications. Wish list? edit: The nice thing is we can produce a lower triangular matrix directly!! from Cholesky based on uBLAS. Returning a full matrix is not wrong but the code is not optimal.

Cuchulainn · December 1st, 2011, 6:20 am

I had a look at some libraries (e.g. Numerical Recipes, QL) and It seems that the following are used:1. Direct matrix solivers (LU, Thomas algo, Gaussian elimination).2. Iterative solvers (Gauss-Seidel, Jacobi, SOR, SSOR, PSOR, CGM).3. Least squares4. Eigensystems (Gerschgorin, Power,..)5. Unitary tranformation, Householder etc.6. SVD.7. Matrix error utilities (e.g. for postprocessing).(8. Sparse matrix utility). To get it working, I suggest writing these with uBLAS as substrate (uBLAS has a lot of useful plumbing s/w). It is possible to make it independent of uBLAS which is fine, but method accuracy is the issue now.Question:1. Is list [1,7] enough? Do we need more features?2. Parallel versions(?) e.g.Boost Thread, OpenMP.

Polter · December 1st, 2011, 4:56 pm

Cuch, regarding #1: you might also consider the other linear algebra decompositions available in Eigen (I posted a link few days ago) plus docs comparing the requirements(applicability)/accuracy/speed trade-offs.On the other hand, I'd wait with parallel versions as far as matrix algebra is concerned (similarly to Eigen, we can note that more coarse-grained parallelism, at the client-code level, might be more desirable), just as long as the design doesn't prevent parallelization in the future (we might want to consider data-oriented design, perhaps I should elaborate on that in another thread).

Cuchulainn · December 1st, 2011, 7:55 pm

QuoteCuch, regarding #1: you might also consider the other linear algebra decompositions available in Eigen (I posted a link few days ago) plus docs comparing the requirements(applicability)/accuracy/speed trade-offs.Had a look. The functionality is similar to the list I drew up. Was kind of expecting a bigger lib on the lines of NAG or IMSL. Just a remark.QuoteOn the other hand, I'd wait with parallel versions as far as matrix algebra is concerned (similarly to Eigen, we can note that more coarse-grained parallelism, at the client-code level, might be more desirable), just as long as the design doesn't prevent parallelization in the future (we might want to consider data-oriented design, perhaps I should elaborate on that in another thread). I agree. Fine-grained is not optimal and less portable.

Cuchulainn · December 2nd, 2011, 6:42 am

Quartz, outrunHere is the modified Gram Schmidt method + test cases and it could be used for a more generic solution. It is about 5 lines of code. So it is a good test case.

Polter · December 2nd, 2011, 5:41 pm

QuoteOriginally posted by: CuchulainnHad a look. The functionality is similar to the list I drew up. Was kind of expecting a bigger lib on the lines of NAG or IMSL. Just a remark.Oh, I must have misunderstood -- I was under impression you were interested in prioritization. Sure, if would like to contribute "a bigger lib on the lines of NAG or IMSL" and are fine with it being a priority in entirety, I think that's a great idea, no need to limit yourself!

Cuchulainn · December 3rd, 2011, 7:51 am

QuoteOriginally posted by: PolterQuoteOriginally posted by: CuchulainnHad a look. The functionality is similar to the list I drew up. Was kind of expecting a bigger lib on the lines of NAG or IMSL. Just a remark.Oh, I must have misunderstood -- I was under impression you were interested in prioritization. Sure, if would like to contribute "a bigger lib on the lines of NAG or IMSL" and are fine with it being a priority in entirety, I think that's a great idea, no need to limit yourself! :-)I think the functionality in Eigen satisfies most of the needs quants need? It is similar to what QL has as well. As a rough estimate, it would take about [100, 140] manhours to get the basic functionality operational, especially if we build on top of the extensive functionality in uBLAS.

Cuchulainn · December 11th, 2011, 4:03 pm

What about some standard C++ for Numerical Quadrature?what about this idea?

Cuchulainn · December 13th, 2011, 5:18 am

QuoteOriginally posted by: outrunQuoteOriginally posted by: outrunShall we do that?..maybe we need to change culture a bit!We'll probably move it in the sandbox on the first day we actually run into the need of Quadrature, so I'm parking this discussion (for myself) until it becomes urgent and actual requirement.I had a look at Boost Quad a few years ago and it looked comprehensive (it is based on QUADPACK?). It would be nice to use just like the other libraries.

Cuchulainn · December 18th, 2011, 2:40 pm

QuoteOriginally posted by: outrunQuoteOriginally posted by: outrunShall we do that?..maybe we need to change culture a bit!We'll probably move it in the sandbox on the first day we actually run into the need of Quadrature, so I'm parking this discussion (for myself) until it becomes urgent and actual requirement.It would be nice to kick the tyres. It is just a simple zip file?

Cuchulainn · October 11th, 2012, 2:33 pm

Here is a popular version of nonunform mesher that is used for Heston etc.Please feel free to improve and feedback.Choose various values of alpha. Could we use this with ADE?

FinancialAlex · October 14th, 2012, 4:23 pm

That Mesher is a very good choice, but it clusters the grid only around a single point. Its extension to clustering around multiple points requires numerically solving ODEs (for single point case we can solve analytically that ODE).