element-wise operations with boost c++ ublas matrix and vector types

Question

i'd like to perform element-wise functions on boost matrix and vector types, e.g. take the logarithm of each element, exponentiate each element, apply special functions, such as gamma and digamma, etc. (similar to matlab's treatment of these functions applied to matrices and vectors.)

i suppose writing a helper function that brute-forced this for each desired function would suffice, but this seems wasteful.

likewise, the boost wiki offers some code to vectorize standard functions, but this seems quite complex.

valarray has been suggested, but i'd like to avoid converting between data types, as i need the ublas data types for other operations (matrix products, sparse matrices, etc.)

any help is greatly appreciated.

Answer 1

The use of begin1() / end1() won't work because it provides access to the element in the default column position (0): consequently, you just access all the elements in the first column. It is better (in the sense that you get the expected behavior) to get sequential access via:

std::transform(mat.data().begin(), mat.data().end(),
               mat.data().begin(), boost::math::tgamma) ;

I suspect this may be a case where the implementation is not quite complete.

Enjoy!

Answer 2

WARNING

The following answer is incorrect. See Edit at the bottom. I've left the original answer as is to give context and credit to those who pointed out the error.

---

I'm not particularly familiar with the boost libraries, so there may be a more standard way to do this, but I think you can do what you want with iterators and the STL transform function template. The introduction to the uBLAS library documentation says its classes are designed to be compatible with the same iterator behavior that is used in the STL. The boost matrix and vector templates all have iterators which can be used to access the individual elements. The vector has begin() and end(), and the matrix has begin1(), end1(), begin2(), and end2(). The 1 varieties are column-wise iterators and the 2 varieties are row-wise iterators. See the boost documentation on VectorExpression and MatrixExpression for a little more info.

Using the STL transform algorithm, you can apply a function to each element of an iterable sequence and assign the result to a different iterable sequence of the same length, or the same sequence. So to use this on a boost uBLAS vector you could do this:

using namespace boost::numeric::ublas;

// Create a 30 element vector of doubles
vector<double> vec(30);

// Assign 8.0 to each element.
std::fill(vec.begin(), vec.end(), 8.0);

// Perform the "Gamma" function on each element and assign the result back
// to the original element in the vector.
std::transform(vec.begin(), vec.end(), vec.begin(), boost::math::tgamma);

For a matrix it would be basically the same thing, you would use either the 1 or 2 family of iterators. Which one you choose to use depends on whether the memory layout of your matrix is row major or column major. A cursory scan of the uBLAS documentation leads me to believe that it could be either one, so you will need to examine the code and determine which one is being used so you choose the most efficient iteration order.

matrix<double> mat(30, 30);
.
.
.

std::transform(mat.begin1(), mat.end1(), mat.begin1(), boost::math::tgamma);

The function you pass as the last argument can be a function taking a single double argument and returning a double value. It can also be a functor.

This is not exactly the same as the vectorization example you cited, but it seems like it should be pretty close to what you want.

---

EDIT

Looks like I should have tested my recommendation before making it. As has been pointed out by others, the '1' and '2' iterators only iterate along a single row / column of the matrix. The overview documentation in Boost is seriously misleading on this. It claims that begin1() "Returns a iterator1 pointing to the beginning of the matrix" and that end1() "Returns a iterator1 pointing to the end of the matrix". Would it have killed them to say "a column of the matrix" instead of "matrix"? I assumed that an iterator1 was a column-wise iterator that would iterate over the whole matrix. For the correct way to do this, see Lantern Rouge's answer.