Theory and implementation for the global optimization of a wide class of algorithms is presented via convex/affine relaxations. The basis for the proposed relaxations is the systematic construction of subgradients for the convex relaxations of factorable functions by McCormick [Math. Prog., 10 (1976), pp. 147–175]. Similar to the convex relaxation, the subgradient propagation relies on the recursive application of a few rules, namely, the calculation of subgradients for addition, multiplication, and composition operations. Subgradients at interior points can be calculated for any factorable function for which a McCormick relaxation exists, provided that subgradients are known for the relaxations of the univariate intrinsic functions. For boundary points, additional assumptions are necessary. An automated implementation based on operator overloading is presented, and the calculation of bounds based on affine relaxation is demonstrated for illustrative examples. Two numerical examples for the global optimization of algorithms are presented. In both examples a parameter estimation problem with embedded differential equations is considered. The solution of the differential equations is approximated by algorithms with a fixed number of iterations.