Automated error analysis with correlation
AMDG The topic of treating error analysis has come up on this list several times in the past. The main difficulty is that the uncertainties of variables are not always independent. I recently had an idea of how to overcome this weakness. The first assumption is that all the variables that we start with are independent. Then, each variable keeps track of how it depends on each of the initial variables. The simplest implementation is to create a unique integer id for every new variable constructed. The format of a variable in the follwing discussion is (value, [(id1, alternate_value1), (id2, alternate_value2)]) where the alternate_values are the values that we would get if the value were increased by their uncertainties. value_with_uncertainty x(1.2, 0.2); // holds (1.2, [(0, 1.4)]) value_with_uncertainty y(1.2, 0.2); // holds (1.2, [(1, 1.4)]) Now, adding x + x, we add the values to get (2.4, ...). We then see that both arguments are dependent on the input with id 0. So, we add these two alternate values to get (2.4, [(0, 2.8)]). The total uncertainty is then 0.4. If on the other hand, we add x + y, the value will still be 2.4, but only the first element of the list of uncertainties depends on id 0, so we use the value of y and the adjusted value of x yielding (2.4, [(0, 2.6), ...]). The same rule gives us the dependence on id 1, for a final result of (2.4, [(0, 2.6), (1, 2.6)]) Now, the combined uncertainty is sqrt(0.2*0.2 + 0.2*0.2) = 0.28. There are several ways this can be made more accurate. 1) instead of taking x + uncertainty_of_x we can take several sample points inside the range x +- uncertainty_of_x 2) it should be possible to track dependence on groups of variables in addition to dependence on single variables independently. All of this adds overhead of course, but I think that a reliable mechanism for error propagation is often going to be worth that cost. The code is attached. Comments? In Christ, Steven Watanabe // value_and_uncertainty.hpp // // Copyright (c) 2008 // Steven Watanabe // // Distributed under the Boost Software License, Version 10.0 (See // accompanying file LICENSE_1_0.txt or copy at // http://www.boost.org/LICENSE_1_0.txt) #ifndef VALUE_AND_UNCERTAINTY_HPP_INCLUDED #define VALUE_AND_UNCERTAINTY_HPP_INCLUDED #include <vector> #include <utility> #include <numeric> #include <cmath> #include <boost/lambda/lambda.hpp> #include <boost/lambda/bind.hpp> template<class T> class value_and_uncertainty { public: value_and_uncertainty(const T& val, const T& unc) : value_(val) { if(unc != T()) { alternate_values.push_back(std::make_pair(new_id(), val + unc)); } } template<class F> value_and_uncertainty(F f, const value_and_uncertainty& arg0) : value_(f(arg0.value_)) { typename alternate_t::const_iterator begin = arg0.alternate_values.begin(), end = arg0.alternate_values.end(); for(; begin != end; ++begin) { alternate_values.push_back(std::make_pair(begin->first, f(begin->second))); } } template<class F> value_and_uncertainty(F f, const value_and_uncertainty& arg0, const value_and_uncertainty& arg1) : value_(f(arg0.value_, arg1.value_)) { typename alternate_t::const_iterator begin0 = arg0.alternate_values.begin(), end0 = arg0.alternate_values.end(); typename alternate_t::const_iterator begin1 = arg1.alternate_values.begin(), end1 = arg1.alternate_values.end(); while(true) { if(begin0 == end0) { for(; begin1 != end1; ++begin1) { alternate_values.push_back(std::make_pair(begin1->first, f(arg0.value_, begin1->second))); } break; } else if(begin1 == end1) { for(; begin0 != end0; ++begin0) { alternate_values.push_back(std::make_pair(begin0->first, f(begin0->second, arg1.value_))); } break; } else { if(begin0->first < begin1->first) { alternate_values.push_back(std::make_pair(begin0->first, f(begin0->second, arg1.value_))); ++begin0; } else if(begin1->first < begin0->first) { alternate_values.push_back(std::make_pair(begin1->first, f(arg0.value_, begin1->second))); ++begin1; } else { alternate_values.push_back(std::make_pair(begin0->first, f(begin0->second, begin1->second))); ++begin0; ++begin1; } } } } const T& value() const { return(value_); } T uncertainty() const { namespace bll = boost::lambda; using namespace bll; using std::sqrt; return(sqrt(std::accumulate(alternate_values.begin(), alternate_values.end(), T(), _1 + bll::bind(bll::protect(_1 * _1), (bll::bind(&std::pair<id_type, T>::second, _2) - value_))))); } private: T value_; typedef int id_type; static id_type new_id() { static id_type result = 0; return(result++); } typedef std::vector<std::pair<id_type, T> > alternate_t; std::vector<std::pair<id_type, T> > alternate_values; }; template<class T> value_and_uncertainty<T> operator+(const value_and_uncertainty<T>& arg0, const value_and_uncertainty<T>& arg1) { using namespace boost::lambda; return(value_and_uncertainty<T>(_1 + _2, arg0, arg1)); } template<class T> value_and_uncertainty<T> operator*(const value_and_uncertainty<T>& arg0, const value_and_uncertainty<T>& arg1) { using namespace boost::lambda; return(value_and_uncertainty<T>(_1 * _2, arg0, arg1)); } template<class T> std::ostream& operator<<(std::ostream& os, const value_and_uncertainty<T>& arg) { return(os << arg.value() << " +- " << arg.uncertainty()); } #endif
Steven: You seem to be following a partial differential approach, where you record the partial derivative of an expression wrt to each of the variables that compose it. Your particular approach is using a sparse representation of the vector of derivatives wrt to all the variables in your system. Wrt the evaluation of the uncertainty (sqrt (0.2*0.2 + 0.2*0.2) = 0.28 in your example), this just corresponds to a norm for the vector, here the Euclidean norm, with a twist that you're not dividing by the total number of variables (not even by the number of variables that are actually involved in your expression). (Note that with a finite number of variables, all norms are equivalent anyway.) It's not an especially new idea. If you go to second derivatives, you can do the same with the Jacobian matrix of your expression. The discussion then becomes sparse vs. dense representation. Your idea seems to boil down to using a sparse representation of this otherwise well known aproach. Am I missing something? BTW, error propagation can also be handled by boost.interval. Note, I am not talking about tracking dependencies (it's the same kind of thing that makes square(x) have a smaller interval enclosure than x * x). Although in general, for small errors, there isn't much difference and we've used interval analysis exactly for that purpose, assuming all the variables are always independent (you get an guaranteed enclosure either way, just not as tight as it could be if you had know the dependencies). My apologies if I missed another particularity of your approach. -- Hervé Brönnimann hervebronnimann@mac.com On Apr 13, 2008, at 12:51 AM, Steven Watanabe wrote:
There are several ways this can be made more accurate. 1) instead of taking x + uncertainty_of_x we can take several sample points inside the range x +- uncertainty_of_x 2) it should be possible to track dependence on groups of variables in addition to dependence on single variables independently.
All of this adds overhead of course, but I think that a reliable mechanism for error propagation is often going to be worth that cost.
The code is attached.
Comments?
In Christ, Steven Watanabe
AMDG Hervé Brönnimann wrote:
Steven: You seem to be following a partial differential approach, where you record the partial derivative of an expression wrt to each of the variables that compose it. Your particular approach is using a sparse representation of the vector of derivatives wrt to all the variables in your system. Wrt the evaluation of the uncertainty (sqrt (0.2*0.2 + 0.2*0.2) = 0.28 in your example), this just corresponds to a norm for the vector, here the Euclidean norm, with a twist that you're not dividing by the total number of variables (not even by the number of variables that are actually involved in your expression). (Note that with a finite number of variables, all norms are equivalent anyway.) It's not an especially new idea. If you go to second derivatives, you can do the same with the Jacobian matrix of your expression.
It should in theory be possible to generalize this to an arbitrary number of derivatives, specified as a template parameter, right?
The discussion then becomes sparse vs. dense representation. Your idea seems to boil down to using a sparse representation of this otherwise well known aproach. Am I missing something?
So, I've just reinvented yet another wheel... Ok. What are the problems with this approach. For me the most important thing is that I can hide the messy calculations behind a nice API. The sparse representation makes this easy.
BTW, error propagation can also be handled by boost.interval. Note, I am not talking about tracking dependencies (it's the same kind of thing that makes square(x) have a smaller interval enclosure than x * x). Although in general, for small errors, there isn't much difference and we've used interval analysis exactly for that purpose, assuming all the variables are always independent (you get an guaranteed enclosure either way, just not as tight as it could be if you had know the dependencies).
I happen to need dependency tracking. In Christ, Steven Watanabe
Steven, A couple of general comments: 1. I would appreciate having this type of error analysis in boost, so I encourage development in this area 2. I would like to have "pluggable" methods of error tracking. For example, most of the time I want partial derivatives with the dependencies tracked, but sometimes I want to know the absolute error of a calculation. 3. I've used the units library quite a bit since it was accepted and made my own class for handling errors, but a nice integration of these two ideas would be nice to see in boost too. Brian
AMDG Schrom, Brian T wrote:
Steven,
A couple of general comments:
1. I would appreciate having this type of error analysis in boost, so I encourage development in this area 2. I would like to have "pluggable" methods of error tracking. For example, most of the time I want partial derivatives with the dependencies tracked, but sometimes I want to know the absolute error of a calculation.
I definitely agree. The code I posted originally was just for proof of concept.
3. I've used the units library quite a bit since it was accepted and made my own class for handling errors, but a nice integration of these two ideas would be nice to see in boost too.
Hopefully they should work together without extra effort... In Christ, Steven Watanabe
Le samedi 12 avril 2008 à 21:51 -0700, Steven Watanabe a écrit :
The simplest implementation is to create a unique integer id for every new variable constructed. The format of a variable in the follwing discussion is (value, [(id1, alternate_value1), (id2, alternate_value2)]) where the alternate_values are the values that we would get if the value were increased by their uncertainties.
This approach to reliable computing is called "affine arithmetic". You may also want to look at other approaches called "Taylor arithmetic" (at first-order) or "slope arithmetic", as they are quite close. Obviously, since they are all first-order arithmetic, they all suffer from second-order terms. Yours is no exception. More precisely, what to do when one performs (x + ex) * (y + ey)? You decided to just discard ex*ey and hope it does not matter. There are two other solutions. First, create a new variable z that always corresponds to x*y. That way, you associate ex*ey to z and no information is lost. The drawback is that you will get a combinatorial explosion on the number of variables. Second, have a value that corresponds to all the higher-terms. The usual approach is to rely on an interval there. That way, you only have one additional term. But the drawback is that you can no longer detect dependencies on second-order terms, which tend to matter in iterative processes. The interval may be stored as lower/upper bounds, or as center/radius values, or as magnitude only (since the center is usually close to zero). An intermediate ground between these two solutions is to have second-order terms as new variables (so quadratic number only) while having an interval for dumping all the higher-order terms. A commonly-encountered variant is to add one single new term for each operation, so that first-order error compensations are properly handled. Anyway, notice that, when using an interval as a wastebin for terms, you can also use it to dump the rounding errors encountered during computations. That way, not only your error analysis is "reliable", but it is also guaranteed. Best regards, Guillaume
participants (4)
-
Guillaume Melquiond -
Hervé Brönnimann -
Schrom, Brian T -
Steven Watanabe