Error propagation: Algorithm for the multi-variate case

The algorithm for the multi-variate case is similar to the single-variate case described above, but slightly more complicated. Multi-variate expressions have the added complication that the terms in the summation of Eq. 1 have to be evaluated separately for each independent variable. This means that a DS and a table of measurement errors (ME) per independent variable have to be maintained.

Symbols (variables and constants) in fussy are tagged with a number or a list of numbers (the IDs) and a type. Symbols representing normal variables are of the type VAR and have a single unique ID, a random error and a partial derivative (of value 1) associated with them. Symbols representing sub-expressions are of the type PARTIAL_VAR and have a list of IDs and a corresponding list of random errors and partial derivatives associated with them. List of unique IDs and the random errors of the independent variables in the expression on the right-hand side (RHS) of the assignment operator constitute the list of IDs and random errors for the PARTIAL_VAR type symbols. When a symbol (of either type) is pushed on the VMS, the entire list of associated IDs is copied to the ID list of the object on the stack. The corresponding random errors and partial derivatives are also copied in the appropriate locations in the ME table and pushed on the appropriate DS respectively. For example, let the IDs of $x_1$ and $x_2$ be $1$ and $2$ respectively. The result of $x_1 * x_2$ on the top of the VMS will have an ID list of {1,2} retaining the information that the result is statistically dependent on the independent variates $x_1$ and $x_2$. If this result is further used as part of another expression, this information will be used to propagate the chain rule for these variables correctly.

Since any expression is built using basic mathematical operators or built-in functions, for the purpose of proving the correctness of the error propagation algorithm for any arbitrary expression, it is sufficient to prove that the algorithm works for the fundamental mathematical operators and built-in functions. The algorithms for evaluating the partial derivatives involving mathematical operators and the final evaluation of the resulting error is described below as pseudo code (see Appendix A for an example). The algorithms are described using the following pseudo functions:

• push(S): to push the symbol or value S on the VMS.
• push(S,DS[i]): to push the symbol or value S on the $i^{th}$ DS.
• pop(): to pop a symbol from the VMS.
• Top(DS[i]): represents the value on the top of the $i^{th}$ DS.

The DS is indexed by the symbol ID(s) (say, N). When a symbol from one of the symbol tables is pushed on the VMS, a value of $1$ is pushed on DS[N] and the associated measurement error is copied into ME[N]. L and R in the pseudo code represents the symbols on the left-hand-side (LHS) and right-hand side (RHS) of the operator respectively. Two functions $f(x,a)$ and $g(x,b)$ with one common variable ($x$) and one non-common variable each ($a$ and $b$) are used as the LHS and RHS operands in the explanation below.

All the operators described below are binary operators. They all pop two values from the top of the VMS, compute the result by applying the corresponding operator, store the value in a temporary stack object, set its ID to the union of the IDs of the operands and push it on the VMS. These operations are performed by the following pseudo function:

   ComputeResult(L,R,Expr)
     {
         L = pop();     R = pop();
         S = Expr(L,R);
         S.IDList = union(L.IDList,R.IDList);
         push(S);
     }

Expr implements the arithmetic of the mathematical operation on L and R.

Partial derivatives with respect to all the IDs in the ID lists of the operands are at the top of the corresponding DSs. For all the IDs common between the two operands, two values are poped from the corresponding DS, say dxR and dxL. The common IDs represent the variables which are part of both the operands (here $x$). Operations to compute the partial derivatives with respect to these common variables (equivalent of the ${\partial \over \partial x}$ operator) are represented by the pseudo function dCommonVar below. The function CommonExpr implements the arithmetic for the derivative computation and is set to the appropriate function for the various operators. These values are computed for each ID in the set composed of the intersection of the ID lists of the two operands, and pushed on the corresponding DS (here, the ID of $x$). Finally, all IDs common between the two operands are removed from the ID lists of each operand.

   dCommonVar(L,R,CommonExpr) 
     { 
         IDList = intersection(L.ID,R.ID); 
         for ID in IDList 
           { 
              dxR = pop(DS[ID]); 
              dxL = pop(DS[ID]); 
              push(CommonExpr(L,R,dxL,dxR),DS[ID]);
              L.Remove(ID);  R.Remove(ID); 
           } 
     }

The list of IDs of the two operands now has IDs corresponding to the non-common variables only ($a$ and $b$ here). Operations for the partial derivatives of the operands with respect to these variables is represented by the pseudo function dNonCommonVar below. LExpr and RExpr computes the value of these derivatives for the LHS and RHS of the operator (equivalent of computing $\frac{\partial f(x,a)}{\partial a}$ and $\frac{\partial g(x,b)}{\partial b}$) using the values from the top of the appropriate DS.

   dNonCommonVar(L,R,LExpr, RExpr)
     {
         for ID in L.IDList
            Top(DS[ID]) = LExpr(L,R,Top(DS[ID]));

         for ID in R.IDList
            Top(DS[ID]) = RExpr(L,R,Top(DS[ID]));
     }

The multiplication operator

The following code returns the result of the operator (L*R) on the top of the VMS with an ID equal to L.ID $\cup$ R.ID.

   Expr(L,R) {return L*R};  /* Compute the value for 
   ComputeResult(L,R,Expr);    the multipication operator */

Following code returns the value(s) of the partial derivative(s) (here dxR*L + dxL*R) with respect to all the common variables in the two operands, on the top of the appropriate DSs.

   Expr(L,R,dxL,dxR) {return dxR*L + dxL*R;}
   dCommonVar(L,R,Expr);

The following code returns the value of the partial derivatives with respect to the variables not common between the operands. These values are returned on the top of the DS corresponding to the remaining IDs of the two operands.

   LExpr(L,R,dx) {return dx*L};
   RExpr(L,R,dx) {return dx*R};
   dNoncommonVar(L,R,LExpr,RExpr)

The division operator

As before, the result of the division operator is computed as L/R and returned on the VMS using the following code.

   Expr(L,R) {return L/R};  /* Compute the value for 
   ComputeResult(L,R,Expr);    the division operator */

For all IDs common between L and R, the partial derivative is computes as (R*dxL - L*dxR)/(R*R) and returned on the appropriate DSs using the following code.

   Expr(L,R,dxL,dxR) {return (R*dxL - L*dxR)/(R*R)};
   dCommonVar(L,R,Expr);

This is equivalent to computing

$${\partial \over {\partial x}}\left[{f(x,a) \over g(x,b)}\rig... ...l x}}} - f(x,a) {{\partial g(x,b)} \over {\partial x}}\right]$$ (3)

Next, partial derivatives with respect to each of the non-common variables are computed (here with respect to $a$ as ${1 \over g(x,b)} {\partial f(x,a) \over \partial a}$ and with respect to $b$ as $-{f(x,a) \over {g(x,b)^2}}{\partial g(x,b) \over \partial b}$) and returned on the appropriate DSs by the following code.

   LExpr(L,R,dx) {return dx/R};
   RExpr(L,R,dx) {return -L*dx/(R*R)};
   dNoncommonVar(L,R,LExpr,RExpr);

The addition operator

The first set of operations is same as that for other operators except that the value pushed on the VMS is L+R.

   Expr(L,R) {return L+R};  /* Compute the value for 
   ComputeResult(L,R,Expr);    the addition operator */

Since the partial derivatives of the expression with respect to the non-common variables is already on the appropriate DS, no separate operation is required for these variables. The partial derivatives with respect to the common set of variables is computed as ${\partial f(x,a) \over {\partial x}}$ and ${\partial g(x,b) \over {\partial x}}$. The pseudo code for these operations is:

   Expr(L,R,dxL,dxR) {return (dxL + dxR)};
   dCommonVar(L,R,Expr);

The subtraction operator

The pseudo code for the subtraction operation is functionally same as that for the addition operator.

   Expr(L,R) {return L-R};  /* Compute the value for the 
   ComputeResult(L,R,Expr);    subtraction operator     */

   /* Compute the derivative w.r.t. common variables */
   Expr(L,R,dxL,dxR) {return (dxL - dxR)};
   dCommonVar(Expr,L,R);

This is equivalent of computing ${\partial f(x,a) \over {\partial x}}-{\partial g(x,b) \over {\partial x}}$. In addition to the above operations, the partial derivatives of the RHS operand with respect to all the non-common variables needs to be negated (here $\partial g(x,b) / \partial b$).

   IDList = union(L.IDList, R.IDList);
   for ID in IDList /* unique IDs in L and R */
     if (ID in R.IDList)
       {
         dxR = pop(DS[ID]);
         push(-dxR, DS[ID]);
       }

The power operator

Again, the result of the to-the-power operator ($L^R$), with an ID equal to L.ID $\cup$ R.ID is pushed on the VMS using the code:

   Expr(L,R) {return L^R};  /* Compute the value for 
   ComputeResult(L,R,Expr);    the power operator */

The partial derivative of the expression with respect to the common variable (here $x$) is computed as:

$${\partial \over {\partial x}}\left(f(x,a)^{g(x,b)}\right)=f(... ...x} + \log(f(x,b)) {\partial g(x,b) \over \partial x} \right]$$ (4)

The pseudo code for this operation is:

   Expr(L,R,dxL,dxR) {return (L^R)*((R/L)*dxL + log(L)*dxR)};
   dCommonVar(L,R,Expr);

The partial derivatives with respect to the non-common set of IDs correspond to $g(x,b)f(x,a)^{[g(x,b)-1]} {\partial f(x,a) \over {\partial a}}$ and $f(x,a)^{g(x,b)} \log(f(x,b)) {\partial g(x,b) \over \partial b}$. Partial derivatives of $f$ and $g$ with respect to $a$ and $b$ are computed using the functions LExpr and RExpr. The computed partial derivatives are pushed back on the appropriate DSs. The pseudo code for this operation is:

   LExpr(L,R,dx) {return R*(L^(R-1))*dx};
   RExpr(L,R,dx) {return (L^R)*log(L)*dx};
   dNoncommonVar(L,R,LExpr,RExpr);

At the terminal operators (e.g. the assignment operator '='), a single value (the result of the right hand side of the terminal operator) is poped from the VMS. The propagated error is then computed using the values from the top of all DSs corresponding to the IDs in the ID list of the poped value. The values from these DSs are the partial derivatives of the expression with respect to the various independent variables used in the expression on the right hand side ( $\partial f / \partial x_i$ in Eq. 1). The corresponding measurement errors ($\delta x_i$ in Eq. 1) are in the appropriated locations in the ME table. Using these values, Eq. 1 is evaluated. This is the final propagated error in the expression.