From 30cfab2e0530e5788fef9b5da76e9fefe42ea42b Mon Sep 17 00:00:00 2001 From: Davis King Date: Thu, 4 Dec 2008 03:15:49 +0000 Subject: [PATCH] Added a comparison of the dlib::matrix to Matlab and also just generally improved the example as a whole. --HG-- extra : convert_revision : svn%3Afdd8eb12-d10e-0410-9acb-85c331704f74/trunk%402701 --- examples/matrix_ex.cpp | 154 +++++++++++++++++++++++++++++++++++++++-- 1 file changed, 150 insertions(+), 4 deletions(-) diff --git a/examples/matrix_ex.cpp b/examples/matrix_ex.cpp index 78ccc8e7e..07d04de06 100644 --- a/examples/matrix_ex.cpp +++ b/examples/matrix_ex.cpp @@ -68,7 +68,7 @@ int main() 1.2, 7.8}; - // load these matrices up with their data. Note that you can only load a matrix + // Load these matrices up with their data. Note that you can only load a matrix // with a C style array if the matrix is statically dimensioned as the M and y // matrices are. You couldn't do it for x since x = M_data would be ambiguous. // (e.g. should the data be interpreted as a 3x3 matrix or a 9x1 matrix?) @@ -87,6 +87,17 @@ int main() cout << "M*x - y: \n" << M*x - y << endl; + // Also note that we can create run-time sized column or row vectors like so + matrix runtime_sized_column_vector; + matrix runtime_sized_row_vector; + // and then they are sized by saying + runtime_sized_column_vector.set_size(3); + + // Similarly, the x matrix can be resized by calling set_size(num rows, num columns). For example + x.set_size(3,4); // x now has 3 rows and 4 columns. + + + // The elements of a matrix are accessed using the () operator like so cout << M(0,1) << endl; // The above expression prints out the value 7.4. That is, the value of @@ -119,6 +130,123 @@ int main() + + // --------------------------------- Comparison with MATLAB --------------------------------- + // Here I list a set of Matlab commands and their equivalent expressions using the dlib matrix. + + matrix A, B, C, D, E; + matrix Aint; + matrix Blong; + + // MATLAB: A = eye(3) + A = identity_matrix(3); + + // MATLAB: B = ones(3,4) + B = uniform_matrix(3,4, 1); + + // MATLAB: C = 1.4*A + C = 1.4*A; + + // MATLAB: D = A.*C + D = pointwise_multiply(A,C); + + // MATLAB: E = A * B + E = A*B; + + // MATLAB: E = A + B + E = A + C; + + // MATLAB: E = E' + E = trans(E); // Note that if you want a conjugate transpose then you need to say conj(trans(E)) + + // MATLAB: E = B' * B + E = trans(B)*B; + + double var; + // MATLAB: var = A(1,2) + var = A(0,1); // dlib::matrix is 0 indexed rather than starting at 1 like Matlab. + + // MATLAB: C = round(C) + C = round(C); + + // MATLAB: C = floor(C) + C = floor(C); + + // MATLAB: C = ceil(C) + C = ceil(C); + + // MATLAB: C = diag(B) + C = diag(B); + + // MATLAB: B = cast(A, "int32") + Aint = matrix_cast(A); + + // MATLAB: A = B(1,:) + A = rowm(B,0); + + // MATLAB: A = B(:,1) + A = colm(B,0); + + // MATLAB: A = [1:5]' + Blong = range(1,5); + + // MATLAB: A = [1:5] + Blong = trans(range(1,5)); + + // MATLAB: A = [1:2:5] + Blong = trans(range(1,2,5)); + + // MATLAB: A = B([1:3], [1:2]) + A = subm(B, range(0,2), range(0,1)); + // or equivalently + A = subm(B, rectangle(0,0,1,2)); + + + // MATLAB: A = B([1:3], [1:2:4]) + A = subm(B, range(0,2), range(0,2,3)); + + // MATLAB: B(:,:) = 5 + set_subm(B,get_rect(B)) = 5; + // or equivalently + set_all_elements(B,5); + + + // MATLAB: B([1:2],[1,2]) = 7 + set_subm(B,range(0,1), range(0,1)) = 7; + + // MATLAB: B([1:3],[2:3]) = A + set_subm(B,range(0,2), range(1,2)) = A; + + // MATLAB: B(:,1) = 4 + set_colm(B,0) = 4; + + // MATLAB: B(:,1) = B(:,2) + set_colm(B,0) = colm(B,1); + + // MATLAB: B(1,:) = 4 + set_rowm(B,0) = 4; + + // MATLAB: B(1,:) = B(2,:) + set_rowm(B,0) = rowm(B,1); + + // MATLAB: var = det(E' * E) + var = det(trans(E)*E); + + // MATLAB: C = pinv(E) + C = pinv(E); + + // MATLAB: C = inv(E) + C = inv(E); + + // MATLAB: [A,B,C] = svd(E) + svd(E,A,B,C); + + // MATLAB: A = chol(E,'lower') + A = cholesky_decomposition(E); + + // MATLAB: var = min(min(A)) + var = min(A); + // ------------------------- Template Expressions ----------------------------- // Now I will discuss the "template expressions" technique and how it is // used in the matrix object. First consider the following expression: @@ -139,8 +267,8 @@ int main() you have to pay for the cost of constructing a temporary matrix object and allocating its memory. Then you pay the additional cost of copying it over to x. It also gets worse when you have more complex expressions - such as x = y + y + y + M*y which would involve the creation and copying - of 4 temporary matrices. + such as x = round(y + y + y + M*y) which would involve the creation and copying + of 5 temporary matrices. All these inefficiencies are removed by using the template expressions technique. The exact details of how the technique is performed are well @@ -164,7 +292,7 @@ int main() This technique works for expressions of arbitrary complexity. So if you - typed x = y + y + y + M*y it would involve no temporary matrices being + typed x = round(y + y + y + M*y) it would involve no temporary matrices being created at all. Each operator takes and returns only matrix_exp objects. Thus, no computations are performed until the assignment operator requests the values from the matrix_exp it receives as input. @@ -209,11 +337,29 @@ int main() tmp() just evaluates a matrix_exp and returns a real matrix object. So it does the same thing as the above code that uses Mtemp. + Another example of this would be chains of matrix multiplies. For example: + */ + x = M*M*M*M; + // A much faster version of this expression would be + x = tmp(M*M)*tmp(M*M); + /* + Anyway, the point of the above discussion is that you shouldn't multiply complex matrix expressions. You should instead assign the expression to a matrix object and then use that object in the multiply. This will ensure that your multiplies are always fast. + + + Note however, that the following two expressions are not afflicted with the + above problem: */ + double value1 = trans(y)*M*y; + double value2 = trans(y)*M*M*y; + /* + These expressions can be evaluated without using temporaries or + needlessly recalculating things as in the case of the above + examples. + !*/ }