[eigen] Zeros affect sequence of operation

[Brad Bell]

Custom scalar AD types, for example Adolc
https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html
record an operation sequence and can use it at different values for the
independent variables.

It appears that Eigen will short circuit certain operations when a value
is zero. In CppAD there is a notion of a variable and a parameter.
Parameters do not depend on the independent variables.

Because CppAD can have multiple levels of AD, a parameter at one level
may be a variable at another. Therefore, in CppAD, it is necessary to
have the notation of identically equal zero; i.e., this value will
always be zero no matter what. Multiplication of a variable by an
identically zero value results in an identically zero value (which leads
to certain optimizations).

Note that, for the type double, identically zero and zero are the same.

It would be nice if Eigen understood variables and parameters and chose
its operation sequence accordingly.

Brad.

[Gael Guennebaud]

Hi,

do you have some specific exemples to point out ?

gael

Post by Brad Bell
Custom scalar AD types, for example Adolc
https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html
record an operation sequence and can use it at different values for the
independent variables.
It appears that Eigen will short circuit certain operations when a value
is zero. In CppAD there is a notion of a variable and a parameter.
Parameters do not depend on the independent variables.
Because CppAD can have multiple levels of AD, a parameter at one level may
be a variable at another. Therefore, in CppAD, it is necessary to have the
notation of identically equal zero; i.e., this value will always be zero no
matter what. Multiplication of a variable by an identically zero value
results in an identically zero value (which leads to certain optimizations).
Note that, for the type double, identically zero and zero are the same.
It would be nice if Eigen understood variables and parameters and chose
its operation sequence accordingly.
Brad.

[Brad Bell]

Below is an example program demonstrating that Eigen sometimes drops sparse
matrix entries that have value zero due to the values in the operands
and not due just to the sparsity pattern of the operands.


# include <iostream>
# include <Eigen/Sparse>
int main(void)
{ typedef Eigen::SparseMatrix<double, Eigen::ColMajor> sparse_matrix;
typedef Eigen::SparseLU<sparse_matrix> sparse_solver;
typedef sparse_matrix::InnerIterator column_itr;
//
int n = 5;
sparse_matrix H(n, n), x(n, 1), y(n, 1);
for(int i = 0; i < n; i++)
{ H.insert(i, i) = 1.0; // n by n identity matrix
x.insert(i, 0) = double(i); // transpose of (0, 1, ..., n-1)
}
std::cout << "Use Eigen for sparse computation of y = H * x\n";
y = H * x;
std::cout << "Note that y(0,1) = 0 in inclued in the result\n";
for(int j = 0; j < y.outerSize(); j++)
{ for(column_itr itr(y, j); itr; ++itr)
{ int row = itr.row();
int col = itr.col();
double value = itr.value();
std::cout << "y(" << row << "," << col << ")=" << value << "\n";
}
}
sparse_solver solver;
solver.analyzePattern(H);
solver.factorize(H);
std::cout << "Use Eigen for sparse solution of y = H * x\n";
x = solver.solve(y); // solve H * x = y
std::cout << "Note that x(0,1) = 0 in not inclued in the result\n";
for(int j = 0; j < x.outerSize(); j++)
{ for(column_itr itr(x, j); itr; ++itr)
{ int row = itr.row();
int col = itr.col();
double value = itr.value();
std::cout << "x(" << row << "," << col << ")=" << value << "\n";
}
}
return 0;
}

Post by Gael Guennebaud
Hi,
do you have some specific exemples to point out ?
gael
Custom scalar AD types, for example Adolc
https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html
<https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html>
record an operation sequence and can use it at different values for the independent variables.
It appears that Eigen will short circuit certain operations when a value is zero. In CppAD there
is a notion of a variable and a parameter. Parameters do not depend on the independent variables.
Because CppAD can have multiple levels of AD, a parameter at one level may be a variable at
another. Therefore, in CppAD, it is necessary to have the notation of identically equal zero;
i.e., this value will always be zero no matter what. Multiplication of a variable by an
identically zero value results in an identically zero value (which leads to certain optimizations).
Note that, for the type double, identically zero and zero are the same.
It would be nice if Eigen understood variables and parameters and chose its operation sequence
accordingly.
Brad.

[Gael Guennebaud]

Indeed, there is no solver within Eigen that is able to explicitly handle a
sparse right hand side/sparse result. So when you call x = solver.solve(y),
internally, y is converted to a dense vector, and the dense result is
turned into a sparse x by pruning zeros:

VectorXd y_dense = y;
VectorXd x_dense = solver.solve(y_dense);
x = x_dense.sparseView();

In your application, is x really sparse ? Otherwise better use dense
rhs/results. To properly handle sparse rhs, all we need is to implement a
sparse triangular solver with sparse rhs. This requires a
depth-first-search traversal to find all the non-zero entries of the
result. But the result is quite dense, like 20% of non-zero, then this
version will very likely be much slower.

gael

Post by Brad Bell
Below is an example program demonstrating that Eigen sometimes drops sparse
matrix entries that have value zero due to the values in the operands
and not due just to the sparsity pattern of the operands.
# include <iostream>
# include <Eigen/Sparse>
int main(void)
{ typedef Eigen::SparseMatrix<double, Eigen::ColMajor> sparse_matrix;
typedef Eigen::SparseLU<sparse_matrix> sparse_solver;
typedef sparse_matrix::InnerIterator column_itr;
//
int n = 5;
sparse_matrix H(n, n), x(n, 1), y(n, 1);
for(int i = 0; i < n; i++)
{ H.insert(i, i) = 1.0; // n by n identity matrix
x.insert(i, 0) = double(i); // transpose of (0, 1, ..., n-1)
}
std::cout << "Use Eigen for sparse computation of y = H * x\n";
y = H * x;
std::cout << "Note that y(0,1) = 0 in inclued in the result\n";
for(int j = 0; j < y.outerSize(); j++)
{ for(column_itr itr(y, j); itr; ++itr)
{ int row = itr.row();
int col = itr.col();
double value = itr.value();
std::cout << "y(" << row << "," << col << ")=" << value << "\n";
}
}
sparse_solver solver;
solver.analyzePattern(H);
solver.factorize(H);
std::cout << "Use Eigen for sparse solution of y = H * x\n";
x = solver.solve(y); // solve H * x = y
std::cout << "Note that x(0,1) = 0 in not inclued in the result\n";
for(int j = 0; j < x.outerSize(); j++)
{ for(column_itr itr(x, j); itr; ++itr)
{ int row = itr.row();
int col = itr.col();
double value = itr.value();
std::cout << "x(" << row << "," << col << ")=" << value << "\n";
}
}
return 0;
}

Post by Gael Guennebaud
Hi,
do you have some specific exemples to point out ?
gael
Custom scalar AD types, for example Adolc
https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html
<https://eigen.tuxfamily.org/dox/TopicCustomizing_CustomScalar.html>
record an operation sequence and can use it at different values for
the independent variables.
It appears that Eigen will short circuit certain operations when a
value is zero. In CppAD there
is a notion of a variable and a parameter. Parameters do not depend
on the independent variables.
Because CppAD can have multiple levels of AD, a parameter at one
level may be a variable at
another. Therefore, in CppAD, it is necessary to have the notation of
identically equal zero;
i.e., this value will always be zero no matter what. Multiplication
of a variable by an
identically zero value results in an identically zero value (which
leads to certain optimizations).
Note that, for the type double, identically zero and zero are the same.
It would be nice if Eigen understood variables and parameters and
chose its operation sequence
accordingly.
Brad.