andreadelprete / expokit-cpp Goto Github PK

View Code? Open in Web Editor NEW

4.0 4.0 0.0 16.08 MB

CMake 0.50% C++ 11.60% C 81.68% Python 6.21%

matrix-exponential

expokit-cpp's People

Contributors

Stargazers

Watchers

expokit-cpp's Issues

Parallelization of computeIntegral and computeDoubleIntegral

Hi @olimexsmart ,

what do you think about parallelizing the computation of the two integrals? Do you think the overhead coming from the parallelization would be higher than the time gained? Of course this depends on the time taken to compute the matrix exponentials, which in turn depends on size and norm of the matrix A, but do you have some experience on the "price" of parallelization?

computeExpTimesVector inconsistent results?

I wrote a little test program to compare the results of computeExpTimesVector and compute in MatrixExponential.h

The results are dependent from the parameter vec_squarings, which I have no idea what's its purpose. With values 1 and 2 the results of the two functions seems identical, from 3 and above the results of computeExpTimesVector diverge towards inf. This is also the reason why the results of LDSUtility::ComputeXt and its Python original compute_x_T were not the same.

How I need to use vec_squarings?

Generalizing computation of "expm times vector" to "expm times matrix"

The method computeExpTimesVector of the class MatrixExponential provides better performance when one needs to compute the product between expm and a vector. Here I would like to generalize the algorithm so that it could be used even in cases when we need to compute the product between expm and a matrix (assuming this matrix has less columns than rows of course).

Our algorithm for the computation of the exponential of a matrix A actually computes first the exponential of A divided by 2^s, where s is the number of squarings, chosen based on the 1-norm of A. After that, the resulting matrix is multiplied by itself s times, thus getting back the exponential of A. However, if all we need is the product between e^A and a matrix V with c columns, we could replace those s Square Matrix-Matrix Multiplications SMMM with 2^s Matrix-Matrix Multiplications (MMM). More interestingly, we could have a mix of SMMM and MMM, for instance s-m SMMM (with m<s) and 2^m MMM (m here is equivalent to vec_squarings). The question is: what is the optimal choice of m? The answer depends on the size of the matrix. If we assume that the cost for a SMMM of size n is n^3, and for a MMM is c*n^2, then the total cost for s-m SMMM and 2^m MMM is going to be

cost(m) = (s-m) * n^3 + 2^m * c * n^2 = s*n^3 - m*n^3 + 2^m * c * n^2

Removing the first term (which is independent of m) and dividing by n^2 we get:

f(m) = c*2^m - m*n

Then we can choose the value of m that minimizes f(m) by computing its gradient and setting it to zero

f'(m) = c*2^m*ln(2) - n = 0
2^m = n / (c*ln(2))
m = log_2 (n / (c*ln(2)))
m = ln(n / (c*ln(2))) / ln(2)
m = (ln(n) - ln(c*ln(2))) / ln(2)
m = ln(n)/ln(2) - ln(c*ln(2)) / ln(2)
m = (ln(n) - ln(c))/ln(2) - ln(ln(2)) / ln(2)
m = ln(n/c)/ln(2) - ln(ln(2)) / ln(2)
m = 1.4427 * ln(n/c) + 0.529

C++ error from BalancingMethods definition at compilation

Hi,

I am trying to install expokit-cpp (master branch), but compilation fails due to C++ errors in include/BalancingMethods.hpp

skleff@akanagi:~/libs/expokit-cpp/build$ make -j4
/usr/bin/cmake -H/home/skleff/libs/expokit-cpp -B/home/skleff/libs/expokit-cpp/build --check-build-system CMakeFiles/Makefile.cmake 0
/usr/bin/cmake -E cmake_progress_start /home/skleff/libs/expokit-cpp/build/CMakeFiles /home/skleff/libs/expokit-cpp/build/CMakeFiles/progress.marks
make -f CMakeFiles/Makefile2 all
make[1]: Entering directory '/home/skleff/libs/expokit-cpp/build'
make -f src/CMakeFiles/expokit.dir/build.make src/CMakeFiles/expokit.dir/depend
make[2]: Entering directory '/home/skleff/libs/expokit-cpp/build'
cd /home/skleff/libs/expokit-cpp/build && /usr/bin/cmake -E cmake_depends "Unix Makefiles" /home/skleff/libs/expokit-cpp /home/skleff/libs/expokit-cpp/src /home/skleff/libs/expokit-cpp/build /home/skleff/libs/expokit-cpp/build/src /home/skleff/libs/expokit-cpp/build/src/CMakeFiles/expokit.dir/DependInfo.cmake --color=
Scanning dependencies of target expokit
make[2]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
make -f src/CMakeFiles/expokit.dir/build.make src/CMakeFiles/expokit.dir/build
make[2]: Entering directory '/home/skleff/libs/expokit-cpp/build'
[ 10%] Building C object src/CMakeFiles/expokit.dir/blas.c.o
[ 10%] Building C object src/CMakeFiles/expokit.dir/dgexpv.c.o
[ 10%] Building C object src/CMakeFiles/expokit.dir/dgchbv.c.o
[ 10%] Building C object src/CMakeFiles/expokit.dir/dgpadm.c.o
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/cc -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -Dexpokit_EXPORTS -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -O3 -DNDEBUG -fPIC   -o CMakeFiles/expokit.dir/dgchbv.c.o   -c /home/skleff/libs/expokit-cpp/src/dgchbv.c
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/cc -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -Dexpokit_EXPORTS -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -O3 -DNDEBUG -fPIC   -o CMakeFiles/expokit.dir/blas.c.o   -c /home/skleff/libs/expokit-cpp/src/blas.c
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/cc -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -Dexpokit_EXPORTS -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -O3 -DNDEBUG -fPIC   -o CMakeFiles/expokit.dir/dgexpv.c.o   -c /home/skleff/libs/expokit-cpp/src/dgexpv.c
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/cc -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -Dexpokit_EXPORTS -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -O3 -DNDEBUG -fPIC   -o CMakeFiles/expokit.dir/dgpadm.c.o   -c /home/skleff/libs/expokit-cpp/src/dgpadm.c
[ 13%] Building C object src/CMakeFiles/expokit.dir/clock.c.o
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/cc -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -Dexpokit_EXPORTS -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -O3 -DNDEBUG -fPIC   -o CMakeFiles/expokit.dir/clock.c.o   -c /home/skleff/libs/expokit-cpp/src/clock.c
[ 16%] Building C object src/CMakeFiles/expokit.dir/dgphiv.c.o
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/cc -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -Dexpokit_EXPORTS -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -O3 -DNDEBUG -fPIC   -o CMakeFiles/expokit.dir/dgphiv.c.o   -c /home/skleff/libs/expokit-cpp/src/dgphiv.c
[ 18%] Building C object src/CMakeFiles/expokit.dir/lapack.c.o
[ 21%] Building C object src/CMakeFiles/expokit.dir/mataid.c.o
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/cc -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -Dexpokit_EXPORTS -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -O3 -DNDEBUG -fPIC   -o CMakeFiles/expokit.dir/lapack.c.o   -c /home/skleff/libs/expokit-cpp/src/lapack.c
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/cc -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -Dexpokit_EXPORTS -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -O3 -DNDEBUG -fPIC   -o CMakeFiles/expokit.dir/mataid.c.o   -c /home/skleff/libs/expokit-cpp/src/mataid.c
[ 24%] Building CXX object src/CMakeFiles/expokit.dir/stop-watch.cpp.o
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/c++  -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -Dexpokit_EXPORTS -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -pedantic -Wno-long-long -Wall -Wextra -Wcast-align -Wcast-qual -Wformat -Wwrite-strings -Wconversion  -O3 -DNDEBUG -fPIC   -o CMakeFiles/expokit.dir/stop-watch.cpp.o -c /home/skleff/libs/expokit-cpp/src/stop-watch.cpp
[ 27%] Building CXX object src/CMakeFiles/expokit.dir/statistics.cpp.o
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/c++  -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -Dexpokit_EXPORTS -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -pedantic -Wno-long-long -Wall -Wextra -Wcast-align -Wcast-qual -Wformat -Wwrite-strings -Wconversion  -O3 -DNDEBUG -fPIC   -o CMakeFiles/expokit.dir/statistics.cpp.o -c /home/skleff/libs/expokit-cpp/src/statistics.cpp
[ 29%] Linking CXX shared library libexpokit.so
cd /home/skleff/libs/expokit-cpp/build/src && /usr/bin/cmake -E cmake_link_script CMakeFiles/expokit.dir/link.txt --verbose=1
/usr/bin/c++ -fPIC  -pedantic -Wno-long-long -Wall -Wextra -Wcast-align -Wcast-qual -Wformat -Wwrite-strings -Wconversion  -O3 -DNDEBUG  -L/usr/lib/x86_64-linux-gnu/lib -Wl,-rpath,/usr/lib/x86_64-linux-gnu/lib -lf2c -lm -shared -Wl,-soname,libexpokit.so -o libexpokit.so CMakeFiles/expokit.dir/blas.c.o CMakeFiles/expokit.dir/dgchbv.c.o CMakeFiles/expokit.dir/dgexpv.c.o CMakeFiles/expokit.dir/dgpadm.c.o CMakeFiles/expokit.dir/clock.c.o CMakeFiles/expokit.dir/dgphiv.c.o CMakeFiles/expokit.dir/lapack.c.o CMakeFiles/expokit.dir/mataid.c.o CMakeFiles/expokit.dir/stop-watch.cpp.o CMakeFiles/expokit.dir/statistics.cpp.o -L/usr/lib/x86_64-linux-gnu/lib -Wl,-rpath,/usr/lib/x86_64-linux-gnu/lib -lf2c -lm -Wl,-rpath,/usr/lib/x86_64-linux-gnu/lib 
make[2]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
[ 29%] Built target expokit
make -f tests/CMakeFiles/testLDS.dir/build.make tests/CMakeFiles/testLDS.dir/depend
make -f tests/CMakeFiles/testNoMalloc.dir/build.make tests/CMakeFiles/testNoMalloc.dir/depend
make -f tests/CMakeFiles/testMatrixExp.dir/build.make tests/CMakeFiles/testMatrixExp.dir/depend
make -f tests/CMakeFiles/testBalancing.dir/build.make tests/CMakeFiles/testBalancing.dir/depend
make[2]: Entering directory '/home/skleff/libs/expokit-cpp/build'
cd /home/skleff/libs/expokit-cpp/build && /usr/bin/cmake -E cmake_depends "Unix Makefiles" /home/skleff/libs/expokit-cpp /home/skleff/libs/expokit-cpp/tests /home/skleff/libs/expokit-cpp/build /home/skleff/libs/expokit-cpp/build/tests /home/skleff/libs/expokit-cpp/build/tests/CMakeFiles/testLDS.dir/DependInfo.cmake --color=
make[2]: Entering directory '/home/skleff/libs/expokit-cpp/build'
cd /home/skleff/libs/expokit-cpp/build && /usr/bin/cmake -E cmake_depends "Unix Makefiles" /home/skleff/libs/expokit-cpp /home/skleff/libs/expokit-cpp/tests /home/skleff/libs/expokit-cpp/build /home/skleff/libs/expokit-cpp/build/tests /home/skleff/libs/expokit-cpp/build/tests/CMakeFiles/testNoMalloc.dir/DependInfo.cmake --color=
make[2]: Entering directory '/home/skleff/libs/expokit-cpp/build'
cd /home/skleff/libs/expokit-cpp/build && /usr/bin/cmake -E cmake_depends "Unix Makefiles" /home/skleff/libs/expokit-cpp /home/skleff/libs/expokit-cpp/tests /home/skleff/libs/expokit-cpp/build /home/skleff/libs/expokit-cpp/build/tests /home/skleff/libs/expokit-cpp/build/tests/CMakeFiles/testMatrixExp.dir/DependInfo.cmake --color=
make[2]: Entering directory '/home/skleff/libs/expokit-cpp/build'
cd /home/skleff/libs/expokit-cpp/build && /usr/bin/cmake -E cmake_depends "Unix Makefiles" /home/skleff/libs/expokit-cpp /home/skleff/libs/expokit-cpp/tests /home/skleff/libs/expokit-cpp/build /home/skleff/libs/expokit-cpp/build/tests /home/skleff/libs/expokit-cpp/build/tests/CMakeFiles/testBalancing.dir/DependInfo.cmake --color=
Scanning dependencies of target testLDS
Scanning dependencies of target testMatrixExp
Scanning dependencies of target testBalancing
Scanning dependencies of target testNoMalloc
make[2]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
make -f tests/CMakeFiles/testBalancing.dir/build.make tests/CMakeFiles/testBalancing.dir/build
make[2]: Entering directory '/home/skleff/libs/expokit-cpp/build'
make[2]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
make -f tests/CMakeFiles/testMatrixExp.dir/build.make tests/CMakeFiles/testMatrixExp.dir/build
make[2]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
make -f tests/CMakeFiles/testNoMalloc.dir/build.make tests/CMakeFiles/testNoMalloc.dir/build
make[2]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
make -f tests/CMakeFiles/testLDS.dir/build.make tests/CMakeFiles/testLDS.dir/build
make[2]: Entering directory '/home/skleff/libs/expokit-cpp/build'
make[2]: Entering directory '/home/skleff/libs/expokit-cpp/build'
make[2]: Entering directory '/home/skleff/libs/expokit-cpp/build'
[ 32%] Building CXX object tests/CMakeFiles/testBalancing.dir/testBalancing.cpp.o
cd /home/skleff/libs/expokit-cpp/build/tests && /usr/bin/c++  -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -pedantic -Wno-long-long -Wall -Wextra -Wcast-align -Wcast-qual -Wformat -Wwrite-strings -Wconversion  -O3 -DNDEBUG    -DBOOST_TEST_MODULE=testBalancingTest -o CMakeFiles/testBalancing.dir/testBalancing.cpp.o -c /home/skleff/libs/expokit-cpp/tests/testBalancing.cpp
[ 35%] Building CXX object tests/CMakeFiles/testMatrixExp.dir/testMatrixExp.cpp.o
cd /home/skleff/libs/expokit-cpp/build/tests && /usr/bin/c++  -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -pedantic -Wno-long-long -Wall -Wextra -Wcast-align -Wcast-qual -Wformat -Wwrite-strings -Wconversion  -O3 -DNDEBUG    -DBOOST_TEST_MODULE=testMatrixExpTest -o CMakeFiles/testMatrixExp.dir/testMatrixExp.cpp.o -c /home/skleff/libs/expokit-cpp/tests/testMatrixExp.cpp
[ 37%] Building CXX object tests/CMakeFiles/testNoMalloc.dir/testNoMalloc.cpp.o
cd /home/skleff/libs/expokit-cpp/build/tests && /usr/bin/c++  -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -pedantic -Wno-long-long -Wall -Wextra -Wcast-align -Wcast-qual -Wformat -Wwrite-strings -Wconversion  -O3 -DNDEBUG    -DBOOST_TEST_MODULE=testNoMallocTest -o CMakeFiles/testNoMalloc.dir/testNoMalloc.cpp.o -c /home/skleff/libs/expokit-cpp/tests/testNoMalloc.cpp
[ 40%] Building CXX object tests/CMakeFiles/testLDS.dir/testLDS.cpp.o
cd /home/skleff/libs/expokit-cpp/build/tests && /usr/bin/c++  -DEIGEN_RUNTIME_NO_MALLOC -DPROFILERYES -I/home/skleff/libs/expokit-cpp/build -I/home/skleff/libs/expokit-cpp/build/include -I/home/skleff/libs/expokit-cpp/include -isystem /opt/openrobots/include/eigen3  -pedantic -Wno-long-long -Wall -Wextra -Wcast-align -Wcast-qual -Wformat -Wwrite-strings -Wconversion  -O3 -DNDEBUG    -DBOOST_TEST_MODULE=testLDSTest -o CMakeFiles/testLDS.dir/testLDS.cpp.o -c /home/skleff/libs/expokit-cpp/tests/testLDS.cpp
In file included from /home/skleff/libs/expokit-cpp/tests/testBalancing.cpp:4:0:
/home/skleff/libs/expokit-cpp/include/BalancingMethods.hpp:204:108: error: no ‘int expokit::BalancingMethods<T, N>::balanceCombined(expokit::BalancingMethods<T, N>::RefMatrix&, expokit::BalancingMethods<T, N>::RefOutMatrix, expokit::BalancingMethods<T, N>::RefOutMatrix, expokit::BalancingMethods<T, N>::RefOutMatrix)’ member function declared in class ‘expokit::BalancingMethods<T, N>’
 int BalancingMethods<T, N>::balanceCombined(RefMatrix& A, RefOutMatrix B, RefOutMatrix D, RefOutMatrix Dinv) {
                                                                                                            ^
In file included from /home/skleff/libs/expokit-cpp/include/MatrixExponential.hpp:22:0,
                 from /home/skleff/libs/expokit-cpp/tests/testMatrixExp.cpp:4:
/home/skleff/libs/expokit-cpp/include/BalancingMethods.hpp:204:108: error: no ‘int expokit::BalancingMethods<T, N>::balanceCombined(expokit::BalancingMethods<T, N>::RefMatrix&, expokit::BalancingMethods<T, N>::RefOutMatrix, expokit::BalancingMethods<T, N>::RefOutMatrix, expokit::BalancingMethods<T, N>::RefOutMatrix)’ member function declared in class ‘expokit::BalancingMethods<T, N>’
 int BalancingMethods<T, N>::balanceCombined(RefMatrix& A, RefOutMatrix B, RefOutMatrix D, RefOutMatrix Dinv) {
                                                                                                            ^
In file included from /home/skleff/libs/expokit-cpp/include/MatrixExponential.hpp:22:0,
                 from /home/skleff/libs/expokit-cpp/include/LDSUtility.hpp:6,
                 from /home/skleff/libs/expokit-cpp/tests/testNoMalloc.cpp:1:
/home/skleff/libs/expokit-cpp/include/BalancingMethods.hpp:204:108: error: no ‘int expokit::BalancingMethods<T, N>::balanceCombined(expokit::BalancingMethods<T, N>::RefMatrix&, expokit::BalancingMethods<T, N>::RefOutMatrix, expokit::BalancingMethods<T, N>::RefOutMatrix, expokit::BalancingMethods<T, N>::RefOutMatrix)’ member function declared in class ‘expokit::BalancingMethods<T, N>’
 int BalancingMethods<T, N>::balanceCombined(RefMatrix& A, RefOutMatrix B, RefOutMatrix D, RefOutMatrix Dinv) {
                                                                                                            ^
In file included from /home/skleff/libs/expokit-cpp/include/MatrixExponential.hpp:22:0,
                 from /home/skleff/libs/expokit-cpp/include/LDSUtility.hpp:6,
                 from /home/skleff/libs/expokit-cpp/tests/testLDS.cpp:4:
/home/skleff/libs/expokit-cpp/include/BalancingMethods.hpp:204:108: error: no ‘int expokit::BalancingMethods<T, N>::balanceCombined(expokit::BalancingMethods<T, N>::RefMatrix&, expokit::BalancingMethods<T, N>::RefOutMatrix, expokit::BalancingMethods<T, N>::RefOutMatrix, expokit::BalancingMethods<T, N>::RefOutMatrix)’ member function declared in class ‘expokit::BalancingMethods<T, N>’
 int BalancingMethods<T, N>::balanceCombined(RefMatrix& A, RefOutMatrix B, RefOutMatrix D, RefOutMatrix Dinv) {
                                                                                                            ^
tests/CMakeFiles/testBalancing.dir/build.make:65: recipe for target 'tests/CMakeFiles/testBalancing.dir/testBalancing.cpp.o' failed
make[2]: *** [tests/CMakeFiles/testBalancing.dir/testBalancing.cpp.o] Error 1
make[2]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
CMakeFiles/Makefile2:1642: recipe for target 'tests/CMakeFiles/testBalancing.dir/all' failed
make[1]: *** [tests/CMakeFiles/testBalancing.dir/all] Error 2
make[1]: *** Waiting for unfinished jobs....
tests/CMakeFiles/testNoMalloc.dir/build.make:65: recipe for target 'tests/CMakeFiles/testNoMalloc.dir/testNoMalloc.cpp.o' failed
make[2]: *** [tests/CMakeFiles/testNoMalloc.dir/testNoMalloc.cpp.o] Error 1
make[2]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
CMakeFiles/Makefile2:1568: recipe for target 'tests/CMakeFiles/testNoMalloc.dir/all' failed
make[1]: *** [tests/CMakeFiles/testNoMalloc.dir/all] Error 2
tests/CMakeFiles/testMatrixExp.dir/build.make:65: recipe for target 'tests/CMakeFiles/testMatrixExp.dir/testMatrixExp.cpp.o' failed
make[2]: *** [tests/CMakeFiles/testMatrixExp.dir/testMatrixExp.cpp.o] Error 1
make[2]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
CMakeFiles/Makefile2:1605: recipe for target 'tests/CMakeFiles/testMatrixExp.dir/all' failed
make[1]: *** [tests/CMakeFiles/testMatrixExp.dir/all] Error 2
tests/CMakeFiles/testLDS.dir/build.make:65: recipe for target 'tests/CMakeFiles/testLDS.dir/testLDS.cpp.o' failed
make[2]: *** [tests/CMakeFiles/testLDS.dir/testLDS.cpp.o] Error 1
make[2]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
CMakeFiles/Makefile2:1531: recipe for target 'tests/CMakeFiles/testLDS.dir/all' failed
make[1]: *** [tests/CMakeFiles/testLDS.dir/all] Error 2
make[1]: Leaving directory '/home/skleff/libs/expokit-cpp/build'
Makefile:143: recipe for target 'all' failed
make: *** [all] Error 2

Adding the defnition of balanceCombined in the .hpp didn't really help since other C++ errors then popped up , this time coming from balanceRodney and balanceNew . Am I doing something wrong here ?

Warm start 2.0

The previous warm-starting algorithm hasn't given good results, so I am investigating a new way for warm-starting the matrix exponential (all very preliminary).

The key idea is that you have already computed e^A and now you have to compute e^(A+B), where |B| << |A| (the norm is the 1 norm). The scaling-and-squaring algorithm relies on Padè approximants, which approximate the exponential as the ratio between two polynomials. Given that the matrices A we are dealing with have a large norm (greater than 5.4), the used polynomials are always of order 13.

e^A = q(A)^-1 * p(A)
p(A) = sum_{i=0}^13 b_i A^i
q(A) = p(-A)

Computing e^(A+B) thus boils down to computing p(A+B) and q(A+B). Can we exploit the previous computations of p(A) and q(A) to do that?

Option 1

We could approximate the high-order (>k) terms of p(A+B) with those of p(A), thus evaluating only the low-order terms of p(A+B). Since |B| << |A| the high-order terms should be well approximated even if we replace A+B with A:

p(A+B) = sum_{i=0}^13 b_i (A+B)^i = sum_{i=0}^k b_i (A+B)^i + sum_{i=k+1}^13 b_i A^i

Option 2

We could use a truncated Taylor expansion of p(A) around A to evaluate p(A+B). For instance we could stop at the first derivative:

p(A+B) = p(A) + p'(A)*B

where the derivative:

p'(A) = sum_{i=1}^13 i * b_i A^(i-1)

can be easily computed in terms of the powers of A, which have already been evaluated for computing p(A). However, in the computation of p(A) not all the powers of A are explicitly computed, but just the power 2, 4 and 6, thanks to an efficient algorithm than evaluates p(A) with only 6 matrix-matrix multiplications (MMM) rather than 13. If we used the same algorithm to evaluate p'(A), exploiting the fact that the first three multiplications have been already computed because they are the ones needed to compute the three above-mentioned powers of A, we can get p'(A) with at most 3 MMM. To that we need to add another MMM to compute p'(A)*B, getting to a total of 4 MMM, rather than the usual 6 MMM we would need to evaluate p(A+B).

This second option seems more promising to me because the error would be O(B^2), therefore it should work well for small |B|. The main issue I see is that the error would be accumulating over the iterations, so we could only use it for a few iterations before needing to do a standard computation. On the bright side, once we have computed p'(A) we can use it several times, as long as the input matrix is not too far from A. Therefore after the first iterations, which costs 4 MMM, the other ones would cost only 1 MMM, making this method extremely advantageous!

Benchmarks branch results

I needed some test code to understand if the changes to the code produced and advantage in time. In particular one of the changes I was dubious about was having function with a return value (exploiting RVO and NRVO) or a in-out parameter:

void compute(const MatrixType& arg, MatrixType& result);
Matrix<type, N, N> compute(const RefMatrix& arg);

Basically nothing, the data structures we move seems too small to impact in significantly the execution. The most reproducible difference was between PT1 and PT4, 4 contacts, 10000 calls:

ComputeDoubleIntegralXt(A, b, xInit, 1, res) took 3113663 us
res = ComputeDoubleIntegralXt(A, b, xInit, 1) took 3087557 us

Which is about a 0.25us of advantage per call for the return value version.
The test call for 10000 times the function in 10 different runs (thus 10e5 calls) to minimize OS noise.
All of this on an Intel N5000 which doesn't even have L3 cache, to underline the fact that on a CPU with more cache and memory bandwidth this is even more insignificant.

@andreadelprete

f2c dependency installation

i tried compiling http://www.netlib.org/f2c/libf2c.zip by reading the instructions in the included readme file,
it seems to compile, however cmake cannot find.
Is it the correct package ? if so how did you get to install it properly ?

Balancing for the 2009 Scaling & Squaring matrix exponential algorithm

I start a new discussion about a balancing algorithm tailored for the 2009 scaling and squaring (S&S) matrix exponential algorithm, to avoid polluting the already existing issue where we discuss balancing for the 2005 algorithm.

Namespace expokit

I've just noticed that expokit does not have a namespace. Nothing urgent, but eventually we should add one.

Getting Up-To-Date speedUps

I'm moving here the email thread with the last email I've sent:

Hi everyone
I've just rebased the master branch with the useful speedUps
@andreadelprete I have one question: can TV be used in the simulator context? Or just TS can be used? I am asking this because I need to explain the couple of changes needed for using the objects.
Luca

Since this part is still a work in progress, we could keep this thread updated on how to use the new stuff in the simulator (@hammoudbilal )

Automatic tuning of vec_squarings

The method computeExpTimesVector of the class MatrixExponential at the moment takes as input the parameter vec_squarings, which is crucial for the performance of the algorithm. If the parameter is not specified by the user, then the algorithm tries to guess it based on the following reasoning.

Our algorithm for the computation of the exponential of a matrix A actually computes first the exponential of A divided by 2^s, where s is the number of squarings, chosen based on the 1-norm of A. After that, the resulting matrix is multiplied by itself s times, thus getting back the exponential of A. However, if all we need is the product between e^A and a vector v, we could replace those s Matrix-Matrix Multiplications MMM with 2^s Matrix-Vector Multiplications (MVM). More interestingly, we could have a mix of MMM and MVM, for instance s-m MMM (with m<s) and 2^m MVM (m here is equivalent to vec_squarings). The question is: what is the optimal choice of m? The answer depends on the size of the matrix. If we assume that the cost for a MMM of size n is n^3, and for a MVM is n^2, then the total cost for s-m MMM and 2^m MVM is going to be

cost(m) = (s-m) * n^3 + 2^m * n^2 = s*n^3 - m*n^3 + 2^m * n^2

Removing the first term (which is independent of m) and dividing by n^2 we get:

f(m) = 2^m - m*n

Then we can choose the value of m that minimizes f(m) by computing its gradient and setting it to zero

f'(m) = 2^m*ln(2) - n = 0
2^m = n / ln(2)
m = log_2 (n / ln(2))
m = ln(n / ln(2)) / ln(2)
m = (ln(n) - ln(ln(2))) / ln(2)
m = ln(n)/ln(2) - ln(ln(2))) / ln(2)
m = 0.693 * ln(n) + 0.529

which is a function that looks like this: vec squarings.pdf. Clearly the value obtained with the above expression needs to be rounded to the closest integer.

However, I've seen that the current implementation in practice chooses a suboptimal value of m (too large). I think there might be two reasons. First, the current implementation is bugged because it does not use the expression above to compute m, so we should just replace the current for loop with

vec_squarings = 0.693 * ln(n) + 0.529

The other reason why even after fixing the bug we could still be overestimating vec_squarings is that we are overestimating the computational complexity of matrix-matrix multiplications (MMM), thinking it is n^3, whereas in practice it is lower.

I've seen that the Strassen algorithm has a complexity of O(n^2.8). Also, looking at the benchmarks of Eigen it seems that the computation time for MMM doesn't grow with the cube of the matrix, which is probably due to several improvements related to vectorization and memory access optimization (some of which are discussed here).

API interface dynamical version of LDSUtility

@hammoudbilal this is a continuation of #2

As for now the class is intended to be instantiated many times as the possible number of contacts. Then a specific instance will be used depending on the number of active contacts the robot has. This is to avoid reshaping at run-time. @andreadelprete pointed out that you may need something more flexible that allows to change the number of contacts at run-time.

I was thinking to have something like LDSUtilityDyn which will duplicate the code of LDSUtility. Same goes with MatrixExponential. Not nice but it is better than messing with inheritance when we are trying to shave us from execution time.
The main changes would be:

Move problem dimension N from template to constructor
Add method that reshapes the internal structures given a new dimension N

Let me know if you have something different in mind
Luca

Computing integral of expm

Hi @olimexsmart

as I mentioned in #12 , we need the computation of the integral of expm for implementing the QP-based approach to deal with friction cones in consim. Computing this matrix is rather simple, and I had already done it in the python script utils_LDS_integral.py. Basically you just need to compute the expm of a matrix that is twice as big as A, and that contains only A, zeros and ones. Then you take the top-left corner of the resulting expm. Looking at the script should be enough to understand.

Before doing anything else, I think we should provide Bilal with this method, otherwise he's blocked and cannot move forward with the development of consim. Could you implement it (in c++) here in expokit and then tell Bilal how to use it?

After this is done, we'll work on making this computation more efficient using the approach described in the pdf I posted in #12 .

Balancing & time scaling

Time scaling is currently giving great results, showing a gain of about 40% in computation time, but the way in which we are tuning the scaling parameter T is not general. It would be nice to design an algorithm to find T automatically and that can work in any case.

BALANCING

Moreover, time scaling corresponds to a special case of balancing, that consists in using a similarity transformation given by a diagonal matrix D, so that the norm of B = D^-1 * A * D is much smaller than the norm of A. Then rather than computing directly e^A, we compute it as:

e^A = D * e^B * D^-1

Time scaling corresponds to setting diagonal elements of D to 1 for the first half and to T for the second half. Clearly this is suboptimal because we are using only one degree of freedom, whereas D has many more parameters.

The problem of balancing for the matrix exponential has been discussed in this mathworks post, coming to the conclusion that he didn't know whether to recommend balancing in general, and that many problems could be solved just by avoiding overscaling using a new version of the scaling-and-squaring algorithm from 2009.

BALANCING FOR THE MATRIX EXPONENTIAL

The balancing algorithm of LAPACK that was used for the matrix exponential preconditioning, was actually originally thought for eigenvalue problems. This algorithm minimizes the Frobenious norm of B, which is not really what we care about when computing matrix exponentials. It'd be better to design a new balancing algorithm that minimizes the 1-norm of B, which is what affects the computation of e^B.

linking error "libf2c.so: undefined reference to `MAIN__'"

When trying to compile the library (master branch) I get the following linking error:

[ 70%] Linking CXX executable testIntegral
cd /home/student/devel/src/expokit/cpp/_build/tests && /usr/bin/cmake -E cmake_link_script CMakeFiles/testIntegral.dir/link.txt --verbose=1
/usr/bin/c++    -pedantic -Wno-long-long -Wall -Wextra -Wcast-align -Wcast-qual -Wformat -Wwrite-strings -Wconversion -std=c++11 -O3 -DNDEBUG   CMakeFiles/testIntegral.dir/testIntegral.cpp.o  -o testIntegral -rdynamic ../src/libexpokit.so -L/usr/lib/x86_64-linux-gnu/lib -Wl,-rpath,/usr/lib/x86_64-linux-gnu/lib -lf2c -lm -Wl,-rpath,/home/student/devel/src/expokit/cpp/_build/src 
/usr/lib/gcc/x86_64-linux-gnu/5/../../../x86_64-linux-gnu/libf2c.so: undefined reference to `MAIN__'
collect2: error: ld returned 1 exit status
tests/CMakeFiles/testIntegral.dir/build.make:98: recipe for target 'tests/testIntegral' failed
make[2]: *** [tests/testIntegral] Error 1

Any idea where this comes from?

2011 COMPUTING THE ACTION OF THE MATRIX EXPONENTIAL

I found another paper by AL-MOHY and HIGHAM which seems very relavant for what we're doing. Luckily they didn't do the same thing that us, but for sure it's a paper we should discuss in our state of the art.

Max number of matrix-matrix multiplications in expm

It could be useful to allow the user to specify a maximum number of matrix-matrix multiplications in the computation of the matrix exponential.

Recommend Projects

React

A declarative, efficient, and flexible JavaScript library for building user interfaces.
Vue.js

🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
Typescript

TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
TensorFlow

An Open Source Machine Learning Framework for Everyone
Django

The Web framework for perfectionists with deadlines.
Laravel

A PHP framework for web artisans
D3

Bring data to life with SVG, Canvas and HTML. 📊📈🎉

Recommend Topics

javascript

JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
web

Some thing interesting about web. New door for the world.
server

A server is a program made to process requests and deliver data to clients.
Machine learning

Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
Visualization

Some thing interesting about visualization, use data art
Game

Some thing interesting about game, make everyone happy.

Recommend Org

Facebook

We are working to build community through open source technology. NB: members must have two-factor auth.
Microsoft

Open source projects and samples from Microsoft.
Google

Google ❤️ Open Source for everyone.
Alibaba

Alibaba Open Source for everyone
D3

Data-Driven Documents codes.
Tencent

China tencent open source team.