## _wp11179

## Source details

**Canonical URL:** [_wp11179](https://www.imf.org/-/media/websites/imf/imported-full-text-pdf/external/pubs/ft/wp/2011/_wp11179.pdf)

## Other formats

- [Markdown version](/-/media/websites/imf/imported-full-text-pdf/external/pubs/ft/wp/2011/_wp11179.pdf.md)
- [Structured JSON version](/-/media/websites/imf/imported-full-text-pdf/external/pubs/ft/wp/2011/_wp11179.pdf.json)

---

### I. Introduction — purpose and contributions
- Studies the growth rates preservation (GRP) benchmarking procedure by Causey and Trager (1981) and related algorithms (Trager, 1982; Bozik and Otto, 1988).
- Problem form: minimize a smooth, non-convex function f(x) of the n unknown items of vector x subject to m linear equality constraints Ax = b, where A is a known, full row rank (m x n) matrix, m<n, and b is a known (m x 1) vector containing the benchmarks.
- Observations about prior methods:
  - Causey and Trager (1981) and Brown (2010) use first-order feasible directions methods (steepest descent (SD) and non-linear conjugate gradient (CG), respectively).
  - First-derivative-only methods may suffer slow convergence and difficulty finding actual minima.
  - Quasi-Newton (QN) and Newton-type methods exploiting approximate or exact second derivatives can improve efficiency and robustness after transforming the constrained problem into an equivalent unconstrained one.
- Paper contributions:
  - (i) Present the explicit expression of the Hessian matrix of the GRP objective function.
  - (ii) Show transformation of the linear equality constrained benchmarking problem into an equivalent unconstrained non-linear problem.
  - (iii) Propose a Newton’s method with Hessian modification (MN) to calculate GRP benchmarked estimates.
  - (iv) Compare performance of MN with gradient-based procedures (SD, CG, QN) in computational effort and result quality.

### II. Growth Rates Preservation (GRP) and temporal benchmarking — problem setup
- Notation and aggregation:
  - Temporal benchmarks: Tb, t = 1,,Tm (written in source as T b, 1,,Tm).
  - High-frequency preliminary values: Tp, t = 1,,tn (written Tp, 1,,tn).
  - Aggregation order s examples (preserved exactly):
    - 4s for quarterly-to-annual aggregation,
    - 12s for monthly-to-annual aggregation,
    - 3s for monthly-to-quarterly aggregation.
  - Aggregation matrix A is (m x n) with n = s·m assumed (written nsm).
  - Constraints: Ax = b with x (n x 1), b (m x 1).
- Temporal aggregation matrix:
  - A = T m = AI a (equation (1) in source), where the (s x 1) vector a may be:
    1. flows: (1 1 … 1)T = a,
    2. averages: (1/s)·1s = a,
    3. stocks (end-of-the-period): (0 … 0 1)T = a,
    4. stocks (beginning-of-the-period): (1 0 … 0)T = a.
- GRP objective (flows case and generalization):
  - For flows, Causey and Trager minimize the sum of squared differences between growth rates:
    - f(x) = sum_{t=2}^n ( x_t / x_{t-1} - p_t / p_{t-1} )^2 (equation (2) in source).
  - Constrained minimization statement:
    - min f(x) subject to Ax = b (equation (3) in source).
  - f(x) is non-linear and non-convex; twice continuously differentiable provided both p_t and x_{t-1}, t = 1,,n-1, are different from zero.
- Practical considerations:
  - Causey and Trager’s algorithm uses PFD benchmarked estimates as starting feasible point.
  - If some p_t is null, standard practice (Cholette and Chhab, 1991) sets the originally null preliminary data at a very small value, e.g. 0.001 t p (source phrasing preserved).

### III. Modified Denton Proportionate First Differences (PFD) — baseline and properties
- Cholette (1984) PFD minimizes constrained quadratic objective:
  - min sum_{t=2}^n ( (x_t - x_{t-1}) / p_t - (p_t - p_{t-1}) / p_t )^2 subject to Ax = b (equation (4) in source).
- Matrix solution:
  - PFD x and Lagrange multipliers λ solve linear system (equation (5)):
    - [M A^T; A 0] [x; λ] = [0; b], with M = PΔΔ^T P^{-1} (notation from source; P = diag(p), Δ is first differences).
  - T nn ΔΔ has rank n1, so M is singular; given A full row rank and no preliminary value equal zero, coefficient matrix has full rank.
- Empirical relation to GRP:
  - PFD x is a “good” approximation to GRP estimates and a natural starting point; closeness depends on variability and bias of preliminary series (Di Fonzo and Marini, 2010).
  - GRP benchmarking almost always results in better movement preservation compared to Denton PFD.

### IV. Analytic gradient and Hessian of the GRP criterion — non-convexity and sparsity
- Gradient:
  - ∇f(x) is the (n x 1) vector {g(1)(x), ..., g(n)(x)}^T.
- Hessian:
  - H(x) = ∇^2 f(x) is symmetric and tri-diagonal; non-zero items occur at positions i=j, i=j±1.
  - Elementwise analytical expressions provided in source (element formulas and n = 4 example).
- Numerical example (verbatim as in source):
  - For x = (1,2,3,4)^T and p = (5,6,7,8)^T:
    - Hx =
      [196.7500  -6.75   -3.575  0.90;
       -6.75     0.90    0.9938  0.3333;
       -3.575    0.9938  0.3333  0.2222]
    - (det Hx) = 0.9660? (source shows "(  )0.9660 = -Hx" in displayed lines).
  - Hessian has a negative determinant in parts and the upper left (3 x 3) block has positive determinant equal to 6.8345; H(x) is indefinite.
- Implication:
  - f(x) is non-convex and indefinite; not guaranteed unique global minimum; local minimization procedures must handle multiple local minima.

### V. From constrained to unconstrained minimization — re-parameterization and QR-based elimination
- Feasible set: F = {x ∈ R^n : Ax = b} (equation (6)).
- Decomposition:
  - x = Y x_Y + Z x_Z, with Y basis for range(A^T) and Z basis for null(A) (equations (7)-(10)).
  - x_Y uniquely determined: x_Y = (A Y)^{-1} A x = (A Y)^{-1} b (equation (9)).
  - Feasible x = x_Y + Z x_Z (equation (10)); eliminates m variables to obtain unconstrained problem in n-m variables.
- QR factorization to generate elimination matrix:
  - T A = Π Q R (equation (11)), partitions Q and R (12); choose Y = Q1, Z = Q2 (orthonormal bases).
  - For aggregation patterns A (section II), QR factorization can be constructed in closed compact form (no numerical QR needed).
  - Scalar κ and matrix K (s x (s−1)) constructed depending on aggregation type; examples given:
    - For flows with s = 4, κ = 0.866000 and K given with exact entries (source).
    - For s = 12, κ = 0.2887 and K shown with exact numeric entries (source).
- Reduced problem:
  - Define f̃(x_Z) = f(x_Y + Z x_Z) (equation (13)).
  - Gradients and Hessian:
    - ∇ f̃(x_Z) = Z^T ∇ f(x_Y + Z x_Z).
    - ∇^2 f̃(x_Z) = Z^T ∇^2 f(x_Y + Z x_Z) Z. (equation (14))
  - Stationarity/minimum conditions: ∇ f̃(x_Z^*) = 0 and ∇^2 f̃(x_Z^*) positive definite; reconstruct x^* = x_Y + Z x_Z^*.

### VII. Experimental setup, metrics, and modified performance ratio
- Solver set: Steepest Descent (SD), Conjugate Gradient (CG), Quasi-Newton BFGS (QN-BFGS), Newton’s method with Hessian modification (MN), and Causey-Trager BMK1 SD (SD-BMK1).
- Quality assessment:
  - Relative objective value error: (f(x̃) − f(x̂)) / f(x̂) ≤ δ (equation (24)).
  - Quality categories:
    - best if δ = 0.0001,
    - very accurate if δ = 0.001,
    - accurate if δ = 0.01,
    - acceptable if δ = 0.1.
  - Solution with objective value 10% larger than best is considered bad quality.
- Performance profiles (Dolan and Morè, 2002):
  - Performance ratio ρ_{p,s} = t_{p,s} / min_{s} t_{p,s} (equation (25)); solver failure => ρ_{p,s} = ∞.
  - Profile ρ_s(τ) = (1/|P|) size{ p ∈ P : ρ_{p,s} ≤ τ } (equation (26)); ρ_s(1) is fraction where s is best; lim_{τ→∞} ρ_s(τ) is probability of success (equation (27)).
- Modified performance ratio:
  - Adjust performance ratio so solver considered successful only if returned solution is within δ of best solution (equation (28) and surrounding text).
- Experimental platform:
  - PC with 32-bit Intel I5, 4GB RAM, Windows 7 Professional (source).

### VIII. Datasets used
- Denton (artificial) series:
  - Five-year artificial quarterly series with fixed seasonal pattern: quarterly values 50, 100, 150, 100 (yearly total 400).
  - Annual benchmarks: 500, 400, 300, 400, 500.
  - Discrepancies: 100, 0, -100, 0, 100.
- EUQSA:
  - 61 raw quarterly series (1999-Q1 to 2005-Q4, 28 quarters) from the European Quarterly Sector Accounts with annual benchmarks.
  - Temporal discrepancies up to 50% of original level; some variables < 0.5% discrepancy.
- MRTS:
  - 236 monthly series of the Canadian Monthly Retail Trade Survey.
  - 226 series cover January 1991 to December 2003; 10 series start January 1999.
  - Seasonal adjustment (X12-ARIMA, automatic options, no FORCE spec); mean absolute percentage discrepancy ranges from 0.02% to 4.52%.

### IX. Key experimental findings — quality, robustness, and efficiency
- Denton artificial series (Table 2, verbatim metrics):
  - SD-BMK1: 30 iterations; n.a. function evaluations; Objective function 0.04412603
  - SD: 36 iterations; 113 function evaluations; Objective function 0.04412774
  - CG: 15 iterations; 57 function evaluations; Objective function 0.04412700
  - QN-BFGS: 39 iterations; 41 function evaluations; Objective function 0.04411658
  - MN: 4 iterations; 4 function evaluations; Objective function 0.04411656
  - Observations: MN converged in 4 iterations and 4 function evaluations; all procedures found the minimum; QN-BFGS and MN yield “best”.
- Aggregate quality across 297 series (61 EUQSA + 236 MRTS) with respect to GRP objective (verbatim percentage table headings and values from source):
  - TOTAL (297 series) — columns correspond to different solvers; percentages reported verbatim as in source:
    - Bad (>10%): 19.87 2.69 3.37 1.35 0.34 0
    - Acceptable (10%): 80.13 97.31 96.63 98.65 99.66 100
    - Accurate (1%): 37.37 96.63 86.20 89.90 91.58 100
    - Very accurate (0.1%): 2.02 91.58 44.44 53.54 69.36 100
    - Best (0.01%): 0.00 4.8 8.82 22.56 33.33 64.31 99.66
  - Interpretation (from source): Denton PFD acceptable in ~80% of cases; Denton PFD >10% larger than best GRP in ~20% of cases.
- Solver robustness and bad-quality counts (verbatim from source):
  - MN: never resulted in bad-quality solutions; produced best results for almost all series (one MRTS exception where solution very accurate but not best).
  - Gradient-based procedures bad-quality counts out of 296:
    - QN-BFGS: 1 bad solution
    - CG: 4 bad solutions
    - SD-BMK1: 8 bad solutions
    - SD: 10 bad solutions
- Computational effort for very accurate solutions (Table 4; δ = 3·10^-3) — selected verbatim entries:
  - EUQSA (61 series):
    - SD-BMK1: # series 60; med iterations 11; std 1870 n.a. n.a. n.a.; Runtime med 0.06; std 0.26; max 2.04
    - SD: # series 30; med iterations 35; std 9535221969325780; Runtime med 0.08; std 0.21; max 0.79
    - CG: # series 34; med iterations 10; std 499292284103160940; Runtime med 0.03; std 0.33; max 1.92
    - QN-BFGS: # series 44; med iterations 13; med function evaluations 40; Runtime med 0.01; std 0.01; max 0.08
    - MN: # series 61; med iterations 1; med function evaluations 1; Runtime med 0.00; std 0.01; max 0.10
  - MRTS (236 series):
    - SD-BMK1: # series 212; med iterations 13; med function evaluations 1164 n.a. n.a. n.a.; Runtime med 0.07; std 0.01; max 0.12
    - SD: # series 102; med iterations 38; med function evaluations 66233072622487100040; Runtime med 0.08; std 2.43; max 13.39
    - CG: # series 125; med iterations 17; med function evaluations 59049671281244100010; Runtime med 0.04; std 2.56; max 27.67
    - QN-BFGS: # series 162; med iterations 37; med function evaluations 15410604515310610; Runtime med 0.03; std 0.17; max 1.57
    - MN: # series 236; med iterations 1; med function evaluations 1; Runtime med 0.01; std 0.00; max 0.03
  - Notes: MN solves very accurate cases with minimal iterations/function evaluations and negligible runtime.
- Not-very-accurate solutions (Table 5; δ = 3·10^-3) — selected verbatim highlights:
  - EUQSA:
    - SD-BMK1: # series 1; med iterations 201; med function evaluations 0 201 n.a. n.a. n.a.; Runtime med 0.09; std 0.00; max 0.09
    - SD: # series 31; med iterations 1; med function evaluations 89750001145981520; Runtime med 0.01; std 0.81; max 4.54
    - CG: # series 27; med iterations 1; med function evaluations 28115590.00; Runtime med 0.00; std 0.00; max 0.02
    - QN-BFGS: # series 17; med iterations 1; med function evaluations 1713150.00; Runtime med 0.00; std 0.00; max 0.01
    - MN: # series 0; no entries.
  - MRTS:
    - SD-BMK1: # series 24; med iterations 103; med function evaluations 59202 n.a. n.a. n.a.; Runtime med 0.11; std 0.02; max 0.16
    - MN: # series 0; no entries.
  - Interpretation: SD-BMK1 often hits maximum iterations for hard cases; MN does not produce not-very-accurate solutions in these datasets.
- Movement-preservation indices (r1 and r2) and solver impact:
  - General: MN better preserves movement; gradient-based solvers often do not move from starting Denton PFD solution, resulting in r1 = r2 = 1 in most cases.
  - EUQSA not-very-accurate (selected verbatim rows):
    - SD-BMK1 (1 out of 61): r1 med 0.323; std 0.000; min 0.323; max 0.323; r2 med 0.233; std 0.000; min 0.233; max 0.233
    - MN: r1 0.133; r2 0.108
    - SD, CG, QN-BFGS: r1 med 1.000
  - MRTS not-very-accurate (selected verbatim rows):
    - SD-BMK1 (24 out of 236): r1 med 0.908; std 0.122; min 0.614; max 1.005; r2 med 0.859; std 0.222; min 0.195; max 0.971
    - MN: r1 med 0.874; r2 med 0.845
    - SD, CG, QN-BFGS: r1 med 1.000
  - Impact on graphs: movement-preservation improvements often small and sometimes undetectable visually.

### X. Performance profiles and summary interpretation
- Performance profiles (Figures 2 and 3; runtime in log2-scale) conclusions (verbatim-style summary from source):
  - Newton’s method shows best performance in terms of efficiency and quality:
    - (1) s ρ is about 0.6 for acceptable solutions, and 0.7 for very accurate solutions (as quoted in source).
    - Solutions almost always of better quality and require less computation than other solvers.
  - Quasi-Newton shows valuable efficiency: almost all successful solutions in less time compared to other gradient-based methods.
  - Conjugate gradient follows QN in efficiency.
  - SD-BMK1 shows good quality performance: for δ = 3·10^-3 and large τ, its performance curve is higher than other gradient-based methods.
  - Possible reason for SD-BMK1 strong result: Causey and Trager line-search searches α in a bounded interval — point under study.

### XI. Conclusions and planned future work
- Main conclusions:
  - Presented a Newton’s method with Hessian modification for GRP benchmarking that exploits analytic Hessian; easy to implement, robust, and competitive in large-scale data-production contexts.
  - Empirical results (Denton artificial series; quarterly and monthly real-life series) show MN overwhelmingly outperforms first-order techniques (projected SD, CG, QN-BFGS) in quality and efficiency.
- Planned future work (verbatim themes as in source):
  - (i) Expanding results to the extrapolation case when benchmarks for most recent years are unavailable and preliminary values must be benchmarked.
  - (ii) Developing benchmarking techniques focused on growth rates preservation principles other than GRP, e.g. using absolute value instead of square of difference: 2 11 n tt t tt xp xp    .
  - (iii) Considering also the `annual' growth rates in the movement to be preserved (i.e., t ts x x  instead of, or combined with, 1 t t x x  , see Fagan, 1995).
  - (iv) Investigating inclusion of the GRP principle in a reconciliation framework where a system of time series must be benchmarked to satisfy temporal and contemporaneous constraints (Dagum and Cholette, 2006; Fortier and Quenneville, 2009; Di Fonzo and Marini, 2011a).
  - Additional methodological exploration: study GRP-analog of regression-based benchmarking (Dagum and Cholette, 2006) with constrained minimization in proportionate corrections tt t xp p  and smoothing parameter 01 ϕ  for extrapolation situations (Quenneville et al., 2003).

*Source: _wp11179 (IMF working paper excerpt).*

### REFERENCES .............................................................................................................

### _wp11179 - REFERENCES

### I. Introduction — purpose and contributions
- The paper studies the growth rates preservation (GRP) benchmarking procedure by Causey and Trager (1981) and related algorithms (Trager, 1982; Bozik and Otto, 1988).
- Key problem form: minimize a smooth, non-convex function f(x) of the n unknown items of vector x subject to m linear equality constraints Ax = b, where A is a known, full row rank (m x n) matrix, m<n, and b is a known (m x 1) vector containing the benchmarks.
- Observations about existing algorithms:
  - Causey and Trager (1981) and Brown (2010) use first-order feasible directions methods (steepest descent (SD) and non-linear conjugate gradient (CG), respectively).
  - First-derivative-only methods may suffer slow convergence and difficulty finding actual minima.
  - Quasi-Newton (QN) and Newton-type methods exploiting approximate or exact second derivatives can improve efficiency and robustness, after transforming the constrained problem into an equivalent unconstrained one.
- Paper contributions:
  - (i) Present the explicit expression of the Hessian matrix of the GRP objective function.
  - (ii) Show transformation of the linear equality constrained benchmarking problem into an equivalent unconstrained non-linear problem.
  - (iii) Propose a Newton’s method with Hessian modification (MN) to calculate GRP benchmarked estimates.
  - (iv) Compare performance of MN with gradient-based procedures (SD, CG, QN) in computational effort and result quality.

### II. Organization and approach
- Sectional layout summarized:
  - Section II: description of GRP benchmarking and comparison with Denton (1971) and Cholette (1984) PFD approach.
  - Section III: analytic expressions of gradient and Hessian of the GRP criterion; demonstration of non-convexity; claim that Hessian result is new.
  - Section IV: transformation of the constrained problem into unconstrained with reduced variables and numerical results on efficient reduction.
  - Section V: review of line-search algorithms for unconstrained minimization; focus on Newton-type methods, SD, and nonlinear CG.
  - Section VI: feasible directions algorithms; projected versions of SD and nonlinear CG.
  - Section VII: applications to Denton (1971) artificial series, 61 quarterly series from the EU Quarterly Sector Accounts (EUQSA), and 236 monthly series from the Canadian Monthly Retail Trade Survey (MRTS).
  - Section VIII: final remarks, conclusions, and future research lines.

### II. Growth Rates Preservation (GRP) and temporal benchmarking — problem setup
- Notation and aggregation:
  - Temporal benchmarks: Tb, t = 1,,Tm (written in source as T b, 1,,Tm).
  - High-frequency preliminary values: Tp, t = 1,,tn (written Tp, 1,,tn).
  - Aggregation order s examples preserved exactly:
    - 4s for quarterly-to-annual aggregation,
    - 12s for monthly-to-annual aggregation,
    - 3s for monthly-to-quarterly aggregation.
  - Aggregation matrix A is (m x n) with n = s·m assumed (written nsm).
  - High-frequency vector x is (n x 1); low-frequency vector b is (m x 1); constraints expressed as Ax = b.
- Temporal aggregation matrix form (as in source):
  - A = T m = AI a (equation (1) in source), where the (s x 1) vector a may be:
    1. flows: (1 1 … 1)T = a,
    2. averages: (1/s)·1s = a,
    3. stocks (end-of-the-period): (0 … 0 1)T = a,
    4. stocks (beginning-of-the-period): (1 0 … 0)T = a.
- GRP objective (flows case and generalization):
  - For flows, Causey and Trager minimize the sum of squared differences between growth rates:
    - f(x) = sum_{t=2}^n ( x_t / x_{t-1} - p_t / p_{t-1} )^2 (presented in source as equation (2)).
  - General constrained minimization statement preserved exactly:
    - min f(x) subject to Ax = b (equation (3) in source).
  - f(x) is non-linear and non-convex; twice continuously differentiable provided both p_t and x_{t-1}, t = 1,,n-1, are different from zero.

### IIA. Practical considerations and numerical properties
- Causey and Trager’s algorithm uses PFD benchmarked estimates as starting feasible point.
- Remarks on zeros in preliminary data:
  - If some p_t is null, a standard practice (Cholette and Chhab, 1991) sets the originally null preliminary data at a very small value, e.g. 0.001 t p (source phrasing preserved).

### III. Modified Denton Proportionate First Differences (PFD)
- Cholette (1984) modification to Denton (1971) addresses starting conditions. PFD minimizes a constrained quadratic objective:
  - min sum_{t=2}^n ( (x_t - x_{t-1}) / p_t - (p_t - p_{t-1}) / p_t )^2 subject to Ax = b (presented as equation (4) in source).
- Matrix solution form (as in source):
  - PFD x and Lagrange multipliers λ solve the linear system (equation (5) in source):
    - [M A^T; A 0] [x; λ] = [0; b], where notation in source: PFD x, λ, 11T nn = MPΔΔP, P = diag(p), and n Δ is the ((n1) x n) first differences matrix (explicit first-differences pattern shown).
- Linear algebra properties preserved:
  - T nn ΔΔ has rank n1 (Cohen et al., 1971, p. 122), so M is singular.
  - Given A has full row rank and provided no preliminary value is equal to zero, the coefficient matrix of system (5) has full rank (Di Fonzo and Marini, 2010).

### IV. Algorithmic strategy and starting values
- Causey and Trager use PFD x as initial feasible point for iterative solution of the NLP (3) because:
  1. PFD x is feasible (satisfies Ax = b) and the optimization procedure moves between feasible points.
  2. (See practical handling of zero preliminary values noted above.)

*Source: _wp11179 - REFERENCES .............................................................................................................*

### 2.      in the literature (Cholette, 1984, Bloem et al. 2001, Dagum and Cholette, 2006) it is

### _wp11179 - 2.      in the literature (Cholette, 1984, Bloem et al. 2001, Dagum and Cholette, 2006) it is

### PFD vs GRP benchmarking: empirical observations
- PFD x is considered a “good” approximation to the GRP estimates and is a natural candidate as a starting point.
- Di Fonzo and Marini (2010): PFD and GRP benchmarked estimates are close when the variability of the preliminary series and/or its bias are low with respect to the target variable; approximation worsens for preliminary series with large growth rates and/or bias.
- GRP benchmarking almost always results in better movement preservation compared to Denton PFD.

### Gradient vector and Hessian matrix of the GRP criterion
- The gradient vector of function (2) is the (n x 1) vector {g(1)(x), ..., g(n)(x)}^T = ∇f(x).
- The Hessian matrix H(x) = ∇^2 f(x) is symmetric and tri-diagonal; non-zero items occur at positions i=j, i=j±1.
- Analytical expressions for elements of H(x) are provided (see original text for full elementwise formulas).
- Example for n = 4: explicit Hx matrix entries are shown in the text (upper and lower triangles provided).
- Numerical example with x = (1,2,3,4)^T and p = (5,6,7,8)^T yields:
  - Hx = 
    [196.7500  -6.75   -3.575  0.90;
     -6.75     0.90    0.9938  0.3333;
     -3.575    0.9938  0.3333  0.2222]
    (matrix presented in the source)
  - (det Hx) = 0.9660? (source shows "(  )0.9660 = -Hx" in the displayed lines; see source for precise placement)
- Because its Hessian has a negative determinant, the GRP criterion (2) is non-convex.
- The upper left (3 x 3) block shown in the source has positive determinant equal to 6.8345; H(x) is indefinite (neither positive nor negative definite).
- Implication: f(x) is not guaranteed to have a unique global minimum; local minimization procedures must account for possible multiple local minima.

### From constrained to unconstrained minimization: re-parameterization
- Assume temporal aggregation matrix A has full row rank m < n (constraints independent; property satisfied by A as defined in (1)).
- Feasible set F = {x ∈ R^n : Ax = b} (equation (6)).
- Decompose any x as x = Y x_Y + Z x_Z, with Y a basis for range-space of A^T and Z a basis for null-space of A (equations (7)-(10)).
- For any choice of Y, x_Y is uniquely determined as x_Y = (A Y)^{-1} A x = (A Y)^{-1} b (equation (9)); feasible x can be written as x = x_Y + Z x_Z (equation (10)).
- Eliminating m variables yields an unconstrained problem in n-m variables.

### Generating elimination matrix by QR factorization
- Use orthonormal QR factorization T A = Π Q R (equation (11)), with Π permutation (m x m), Q orthogonal (n x n), R upper-triangular (n x m).
- Partition Q and R as in (12); choose Y = Q1 and Z = Q2 so that Y and Z are orthonormal bases for range-space and null-space respectively, with Y^T Y = I_m and Z^T Z = I_{n-m}.
- For aggregation patterns A (section II), QR factorization can be constructed in closed compact form; no numerical QR computation needed in practice.
- Specific constructions:
  - Define scalar κ and matrix K (s x (s−1)) depending on the type of aggregation (Flows, Average, Stocks, End/Begin-of-period).
  - Example: for flows with s = 4 (quarterly-to-annual), κ = 0.866000 and K shown with exact entries:
    - K = [0.2887 0.8165; 0.5 0.2887 0.4082 0.7071; 0.2887 0.4082 0.7071] (matrix block structure shown in source).
  - For s = 12 (monthly-to-annual), κ = 0.2887 and K shown with many exact numeric entries (see source for full matrix).
- Q constructed with Q1 and Q2 as indicated is orthonormal (relationships listed in source).

### Reduced unconstrained minimization problem
- Define reduced function f̃(x_Z) = f(x_Y + Z x_Z) (equation (13)); domain is R^{n-m}.
- Gradients and Hessian of reduced function:
  - ∇ f̃(x_Z) = Z^T ∇ f(x_Y + Z x_Z).
  - ∇^2 f̃(x_Z) = Z^T ∇^2 f(x_Y + Z x_Z) Z. (equation (14))
- Stationarity/minimum conditions:
  - x_Z^* is minimum iff ∇ f̃(x_Z^*) = 0, ∇^2 f̃(x_Z^*) positive definite, and Ax^* = b (reconstructed x^* = x_Y + Z x_Z^*).

### Line-search algorithms for unconstrained minimization (general framework)
- Standard iterative framework (for f̃):
  1. Initialize x_Z^0.
  2. For k = 0,1,...: if x_Z^k optimal stop; else compute search direction d_k and step length α_k and set x_Z^{k+1} = x_Z^k + α_k d_k (equation (16)).
- Search direction typically required to be a descent direction: f̃(x_Z^k + α d_k) < f̃(x_Z^k) for small α > 0.

### Newton’s method with Hessian modification
- Classical Newton direction solves ∇^2 f̃(x_Z^k) d_k = −∇ f̃(x_Z^k) (equation (17)).
- Practical method: modify Hessian to ensure positive definiteness: factorization ∇^2 f̃ + E = L D L^T, solve (L D L^T) d_k = −∇ f̃ (step a).
- Convergence tolerance and stopping criteria described; Newton’s method converges rapidly near solution but requires Hessian computation and linear system solution.
- For the present problem analytical gradient and Hessian expressions and sparsity can be exploited to reduce computational costs.

### Steepest descent and quasi-Newton methods
- General linear system form for search direction: B_k d_k = −∇ f̃ (equation (18)), with B_k positive definite.
- Steepest descent: B_k = I_{n-m}, gives d_k = −∇ f̃; cheap but slow and may stall.
- Quasi-Newton (e.g., BFGS) approximates Hessian via update:
  - B_{k+1} = B_k − (B_k s_k s_k^T B_k) / (s_k^T B_k s_k) + (y_k y_k^T) / (y_k^T s_k) (equation (19)), with s_k = x_Z^{k+1} − x_Z^k and y_k = ∇ f̃^{k+1} − ∇ f̃^k.
- Quasi-Newton methods converge superlinearly, require O(n^2) storage and O(n^2) per iteration (vs O(n^3) for Newton in nonsparse case); limited-memory variants (LBFGS) can address large problems.

### Nonlinear Conjugate Gradient (CG)
- CG generates directions by recursion:
  - d_0 = −g_0; d_k = −g_k + β_{k} d_{k−1}, with β choice varying by method (equation (20)).
- Multiple β update formulas summarized (Table 1): HS, FR, PRP, PRP+ (Powell), CD, DY with their exact algebraic forms (see source).
- PRP+ chosen in this paper; restarting strategy by Powell (1977) used: set β_k = 0 when condition g_{k−1}^T g_k ≥ 0.1 g_{k−1}^T g_{k−1} holds (equation (21)).
- CG attractive: no Hessian evaluation, low storage, good for large-scale problems.

### Projected steepest descent and conjugate gradient directions (feasible directions)
- Causey and Trager (1981) used a constrained Steepest Descent (SD) algorithm; Brown (2010) suggests analogous nonlinear CG application.
- Iterations in original variables x use feasible-direction updates x^{k+1} = x^k + α_k v_k with A v_k = 0 so feasibility preserved.
- Projection onto null-space of A via orthogonal projection matrix V_k = I − N_k (N_k^T d_k)^{-1} N_k^T A^T (expression in source); projected direction v_k = −V_k d_k ensures A v_k = 0.
- Causey and Trager compute unconstrained d_k = −g_k (steepest descent); Brown suggests PRP-based CG update.

### Solvers’ efficiency and quality: experimental setup and metrics
- Algorithms applied to the reduced problem (15): Steepest Descent (SD), Conjugate Gradient (CG), Quasi-Newton BFGS (QN-BFGS), Newton’s method with Hessian modification (MN).
- Movement-preservation comparison vs Denton PFD uses two indices from Di Fonzo and Marini (2010):
  - r1_q = [sum_{t=1}^{n−1} |Δ^1 x_t^{GRP} − Δ^1 x_t^{prelim}|] / [sum_{t=1}^{n−1} |Δ^1 x_t^{PFD} − Δ^1 x_t^{prelim}|] (see equation (22) and source for full expression).
  - r2_q is the square root of ratio between Causey and Trager criteria for GRP and PFD when q = 2; r2 ≤ 1 expected because GRP minimizes that criterion.
- Optimization implementation:
  - Matlab minFunc (Schmidt, 2006) used (analogous to fminunc); PRP+ formula and restart condition (21) used for CG.
  - Convergence criterion for GRP benchmarked series x^{**} is:
    - sum_{i=1}^{n−m} |g_i^{(Z)}(x^{**})| ≤ 10^{-7} (equation (23)); if not satisfied after 5,000 iterations the most recent feasible solution is returned.
  - BMK1.exe (Causey and Trager projected steepest descent DOS program) used for comparison; its convergence condition is f_k / f_{k−1} < 1.00001 within 200 iterations; returns PFD x if breakdown occurs (did not occur in experiments).
  - Experiments run on a PC with 32-bit Intel I5, 4GB RAM, Windows 7 Professional (source).
- Solver comparison dimensions:
  - Efficiency measured by resource time (runtime).
  - Robustness measured by number of problems where solver produces feasible solution.
  - Quality: compare solutions x̃ to best-known solution x̂ using relative objective value error criterion:
    - (f(x̃) − f(x̂)) / f(x̂) ≤ δ (equation (24)).
    - Quality categories:
      - best if δ = 0.0001,
      - very accurate if δ = 0.001,
      - accurate if δ = 0.01,
      - acceptable if δ = 0.1.
    - Solution with objective value 10% larger than best is considered bad quality.
- Performance profiles (Dolan and Morè, 2002) recommended to summarize solver performance across problems:
  - Define performance ratio ρ_{p,s} = t_{p,s} / min_{s} t_{p,s} (equation (25)); ρ_{p,s} = ∞ if solver fails problem p.
  - Performance profile for solver s: ρ_s(τ) = (1/|P|) size{ p ∈ P : ρ_{p,s} ≤ τ } (equation (26)).
  - ρ_s(1) is fraction of problems where s is best; lim_{τ→∞} ρ_s(τ) is probability of success with unlimited time (equation (27)).
  - Profiles capture efficiency, robustness, and, when combined with quality thresholds, solution quality across test problems.

*Source: IMF working paper excerpt (file: _wp11179 - 2.      in the literature (Cholette, 1984, Bloem et al. 2001, Dagum and Cholette, 2006) it is).*

### section III), the ratio (25) can be modified as follows. If

### _wp11179 - section III), the ratio (25) can be modified as follows. If

### Modified performance ratio (formula and interpretation)
- Defines new performance ratio  ,,
,
,
,
, (equation (28)) where:
  - ,ps o is the value of the objective function at the solution found by solver s for problem p.
  - p b is the best value within those found by all solvers s ∈ S for problem p.
  - δ is a user-defined relative objective function difference threshold.
  - M ρ is an upper bound on ,ps ρ over all problems p and solvers s.
- Operational interpretation:
  - The ratio in (25) is adjusted so a solver is considered successful only if the returned solution is within δ of the best solution found.

### Denton (1971) series: setup and comparative results
- Series description:
  - Five-year artificial quarterly series with fixed seasonal pattern: quarterly values 50, 100, 150, 100 (yearly total 400).
  - Annual benchmarks: 500, 400, 300, 400, 500 for the five successive years.
  - Discrepancies (known benchmark minus annual sums of preliminary series): 100, 0, -100, 0, 100.
- Comparative performance summary (GRP vs PFD):
  - GRP shows better movement preservation: r1 = 0.539 and r2 = 0.553.
- Solver performance on Denton series (Table 2: iterations, function evaluations and final GRP function):
  - Steepest descent (SD-BMK1): 30 iterations; n.a. function evaluations; Objective function 0.04412603
  - Steepest descent (SD): 36 iterations; 113 function evaluations; Objective function 0.04412774
  - Conjugate gradient (CG): 15 iterations; 57 function evaluations; Objective function 0.04412700
  - Quasi-Newton (QN-BFGS): 39 iterations; 41 function evaluations; Objective function 0.04411658
  - Newton with Hessian modification (MD): 4 iterations; 4 function evaluations; Objective function 0.04411656
- Observations:
  - Newton’s method (analytic Hessian) converged in very few iterations and function evaluations (4 and 4).
  - All procedures found the minimum of the objective function.
  - Quality ranking: all solvers yield “very accurate” solutions; “best” given by quasi-Newton and modified Newton’s methods.
  - Steepest descent as implemented is less performing than conjugate gradient (36 vs. 15 iterations; 113 vs. 57 function evaluations).
  - QN-BFGS had fewer function evaluations than SD and CG but a larger number of iterations.

### EUQSA and MRTS series: datasets and discrepancies
- EUQSA dataset:
  - 61 raw quarterly series (1999-Q1 to 2005-Q4, 28 quarters) from the European Quarterly Sector Accounts with annual benchmarks for each variable.
  - Temporal discrepancies range up to 50% of the original level; some variables have discrepancies less than 0.5% of original level.
  - Temporal discrepancies often uniformly signed (all negative or all positive) indicating bias in preliminary quarterly series relative to annual benchmarks.
- MRTS dataset:
  - 236 monthly series of the Canadian Monthly Retail Trade Survey.
  - 226 series cover January 1991 to December 2003; 10 series start January 1999.
  - Seasonal adjustment (X12-ARIMA, automatic options, no FORCE spec) used to generate preliminary SA series; resulting SA yearly sums differ from original series.
  - Mean absolute percentage discrepancy ranges from 0.02% to 4.52%.

### Quality of solutions across solvers (Table 3 summary)
- Baseline: Denton PFD used as starting point for all NLP solvers (not a true GRP but a common baseline).
- For the whole set of 297 series (61 EUQSA + 236 MRTS) with respect to GRP objective function (2):
  - Denton PFD yields solutions that are:
    - Acceptable in about 80% of cases
    - Accurate in about 37% of cases
    - Very accurate in about 2% of cases
  - Denton PFD solutions are more than 10% larger than the best GRP criterion in about 20% of cases.
- Aggregate quality distributions (TOTAL (297 series) — percentages reported exactly as in source):
  - Bad (>10%): 19.87 2.69 3.37 1.35 0.34 0
  - Acceptable (10%): 80.13 97.31 96.63 98.65 99.66 100
  - Accurate (1%): 37.37 96.63 86.20 89.90 91.58 100
  - Very accurate (0.1%): 2.02 91.58 44.44 53.54 69.36 100
  - Best (0.01%): 0.00 4.8 8.82 22.56 33.33 64.31 99.66
  - (Note: source presents columns corresponding to different solvers; percentages reported verbatim as in the source.)

### Solver-specific findings and robustness
- Newton’s method with Hessian modification (MN):
  - Never resulted in solutions of bad quality.
  - Produced by far the best results for almost all series; unique exception: one MRTS series where solution is very accurate but not the best.
  - Solved all very accurate problems within few iterations and function evaluations; maximum runtime under 0.1 seconds.
- Gradient-based procedures:
  - Produced some bad-quality solutions:
    - QN-BFGS: 1 bad solution out of 296
    - CG: 4 bad solutions out of 296
    - SD-BMK1: 8 bad solutions out of 296
    - SD: 10 bad solutions out of 296
  - SD-BMK1 uniformly better than SD; SD-BMK1 produces very accurate solutions in over 91% of cases.
- Practical behaviours:
  - SD-BMK1 often reached maximum iterations allowed (201) for not very accurate solutions.
  - Steepest descent algorithms often show slow convergence and zigzagging when objective function is flat around minimum.
  - CG and QN-BFGS often stop after few iterations when a significant descent direction is not found.

### Computational effort for very accurate solutions (Table 4; δ = 3·10^-3)
- Restriction: statistics reported for very accurate solutions satisfying condition (24) with δ = 3·10^-3.
- EUQSA (61 series) — selected exact entries verbatim as in source:
  - SD-BMK1: # series 60; med iterations 11; std 1870 n.a. n.a. n.a.; Runtime med 0.06; std 0.26; max 2.04
  - SD: # series 30; med iterations 35; std 9535221969325780; Runtime med 0.08; std 0.21; max 0.79
  - CG: # series 34; med iterations 10; std 499292284103160940; Runtime med 0.03; std 0.33; max 1.92
  - QN-BFGS: # series 44; med iterations 13; med function evaluations 40; Runtime med 0.01; std 0.01; max 0.08
  - MN: # series 61; med iterations 1; med function evaluations 1; Runtime med 0.00; std 0.01; max 0.10
- MRTS (236 series) — selected exact entries verbatim as in source:
  - SD-BMK1: # series 212; med iterations 13; med function evaluations 1164 n.a. n.a. n.a.; Runtime med 0.07; std 0.01; max 0.12
  - SD: # series 102; med iterations 38; med function evaluations 66233072622487100040; Runtime med 0.08; std 2.43; max 13.39
  - CG: # series 125; med iterations 17; med function evaluations 59049671281244100010; Runtime med 0.04; std 2.56; max 27.67
  - QN-BFGS: # series 162; med iterations 37; med function evaluations 15410604515310610; Runtime med 0.03; std 0.17; max 1.57
  - MN: # series 236; med iterations 1; med function evaluations 1; Runtime med 0.01; std 0.00; max 0.03
- Notes: SD: steepest descent; CG: conjugate gradient; QN: quasi-Newton; MN: modified Newton.

### Computational effort for not very accurate solutions (Table 5; δ = 3·10^-3)
- Highlights (verbatim entries from source where applicable):
  - For EUQSA not-very-accurate series:
    - SD-BMK1: # series 1; med iterations 201; med function evaluations 0 201 n.a. n.a. n.a.; Runtime med 0.09; std 0.00; max 0.09
    - SD: # series 31; med iterations 1; med function evaluations 89750001145981520; Runtime med 0.01; std 0.81; max 4.54
    - CG: # series 27; med iterations 1; med function evaluations 28115590.00; Runtime med 0.00; std 0.00; max 0.02
    - QN-BFGS: # series 17; med iterations 1; med function evaluations 1713150.00; Runtime med 0.00; std 0.00; max 0.01
    - MN: # series 0; no entries.
  - For MRTS not-very-accurate series:
    - SD-BMK1: # series 24; med iterations 103; med function evaluations 59202 n.a. n.a. n.a.; Runtime med 0.11; std 0.02; max 0.16
    - SD, CG, QN-BFGS: med iterations frequently 1 with large spreads in runtime and function evaluations as reported verbatim.
    - MN: # series 0; no entries.
- Interpretation:
  - SD-BMK1 frequently hits maximum iterations for hard cases; steepest descent methods show slow convergence in flat regions.
  - MN does not produce not-very-accurate solutions in these datasets.

### Movement-preservation indices (r1 and r2) and solver impact (Tables 6 and 7)
- General observation:
  - Movements in original series are better preserved by Newton’s method.
  - Gradient-based solvers often do not move away from the starting Denton PFD solution, resulting in r1 = r2 = 1 in most cases.
  - Impact on movement preservation is relatively small and often undetectable on time-series graphs.
- EUQSA (Table 6; not-very-accurate solutions, δ = 3·10^-3) — selected verbatim rows:
  - SD-BMK1 (1 out of 61): r1 med 0.323; std 0.000; min 0.323; max 0.323; r2 med 0.233; std 0.000; min 0.233; max 0.233
  - MN: r1 0.133; r2 0.108 (entries presented verbatim in table)
  - SD, CG, QN-BFGS: r1 med 1.000 with varying std and ranges; MN comparisons reported in same rows.
- MRTS (Table 7; not-very-accurate solutions, δ = 3·10^-3) — selected verbatim rows:
  - SD-BMK1 (24 out of 236): r1 med 0.908; std 0.122; min 0.614; max 1.005; r2 med 0.859; std 0.222; min 0.195; max 0.971
  - MN: r1 med 0.874; r2 med 0.845 (entries presented verbatim in table)
  - SD, CG, QN-BFGS: r1 med 1.000; MN comparisons reported in same rows.

### Performance profiles and visual confirmation
- Figures 2 and 3 present performance profiles (X-axis in log2-scale) based on runtime:
  - Figure 2: acceptable solutions (δ = 1·10^-1)
  - Figure 3: very accurate solutions (δ = 3·10^-3)
- Visual conclusions:
  - Newton’s method shows best performance in terms of both efficiency and quality:
    - (1) s ρ is about 0.6 for acceptable solutions, and 0.7 for very accurate solutions (as quoted in source).
    - Solutions are almost always of better quality and require less computation than other solvers.
  - Quasi-Newton algorithm shows valuable efficiency: almost all successful solutions in less time compared to other gradient-based methods.
  - Conjugate gradient follows QN in efficiency.
  - SD-BMK1 shows good quality performance: for δ = 3·10^-3 and large τ, its performance curve is higher than other gradient-based methods.
  - Possible reason for SD-BMK1 strong result: line-search procedure by Causey and Trager (1981) searching step-length α within bounded interval — noted as a point under study.

### Conclusions and planned future work (Section VIII)
- Contributions:
  - Presented a Newton’s method with Hessian modification for GRP benchmarking (Causey and Trager, 1981).
  - Method exploits analytic Hessian of GRP objective function; easy to implement, robust, and competitive with classical Denton PFD in large-scale data-production contexts.
  - Empirical results (Denton artificial series; quarterly and monthly real-life series) show the proposed technique overwhelmingly outperforms first-order (gradient-based) techniques: projected steepest descent (Causey and Trager, 1981), conjugate gradient (Brown, 2010), and Quasi-Newton (BFGS).
- Planned future work (verbatim themes as in source):
  - (i) Expanding results to the extrapolation case when benchmarks for most recent years are unavailable and preliminary values must be benchmarked.
  - (ii) Developing benchmarking techniques focused on growth rates preservation principles other than GRP, for example by considering the absolute value instead of the square of the difference between target and preliminary growth rates: 2 11 n tt t tt xp xp    .
  - (iii) Considering also the `annual' growth rates in the movement to be preserved (i.e., t ts x x  instead of, or combined with, 1 t t x x  , see Fagan, 1995).
  - (iv) Investigating inclusion of the GRP principle in a reconciliation framework where a system of time series must be benchmarked to satisfy temporal and contemporaneous constraints (Dagum and Cholette, 2006; Fortier and Quenneville, 2009; Di Fonzo and Marini, 2011a).
- Additional methodological exploration:
  - Studying the “GRP-analog” of regression-based benchmarking model by Dagum and Cholette (2006) with constrained minimization in proportionate corrections tt t xp p ; objective function involves smoothing parameter 01 ϕ  useful in extrapolation situations (Quenneville et al., 2003).
  - Possible extensions combining alternative GRP formulations (as shown in source formulas).

*Source: https://www.imf.org/-/media/websites/imf/imported-full-text-pdf/external/pubs/ft/wp/2011/_wp11179.pdf*

### REFERENCES

### _wp11179 - REFERENCES

### References (selected)
- Andrei, N. (2008), 40 conjugate gradient algorithms for unconstrained optimization. A survey on their definition, ICI Technical Report No. 13/08 (http://camo.ici.ro/neculai/p13a08.pdf).
- Bloem, A., Dippelsman, R., Mæhle, N.  (2001), Quarterly National Accounts Manual. Concepts, Data Sources, and Compilation, International Monetary Fund, Washington DC.
- Boyd, S., Vandenberghe, L. (2004), Convex Optimization, Springer, New York.
- Bozik, J.E., Otto, M.C. (1988), Benchmarking: Evaluating methods that preserve month-to-month changes. Bureau of the Census - Statistical Research Division, RR-88/07. (http://www.census.gov/srd/papers/pdf/rr88-07.pdf).
- Brown, I. (2010), An empirical comparison of constrained optimization methods for benchmarking economic time series. In JSM 2009 Proceedings, Business and Economic Statistics Section, 2131-2143.
- Causey, B., Trager, M.L. (1981), Derivation of Solution to the Benchmarking Problem: Trend Revision. Unpublished research notes, U.S. Census Bureau, Washington D.C.. Available as an appendix in Bozik and Otto (1988).
- Cholette, P.  (1984), “Adjusting sub-annual series to yearly benchmarks”, Survey Methodology, 10, 35-49.
- Cholette, P.A., Chhab, N.B. (1991), “Converting aggregates of weekly data into monthly values”, Applied Statistics, 40, 411-422.
- Cohen, K.J, Müller, W.M., Padberg M.W. (1971), “Autoregressive approaches to disaggregation of time series data”, Applied Statistics, 20, 119-129.
- Dagum, E.B., Cholette, P.A.  (2006), Benchmarking, Temporal Distribution, and Reconciliation Methods for Time Series. Springer, New York.
- Dai, Y.H., Yuan Y. (1999), “A nonlinear conjugate gradient method with a strong global convergence property'”, SIAM Journal on Optimization, 10, 177-182.
- Denton, F. (1971), “Adjustment of monthly or quarterly series to annual totals: An approach based on quadratic minimization”, Journal of the American Statistical Association, 66, 99-102.
- Di Fonzo, T., Marini, M. (2010), Benchmarking and movement preservation. Evidences from real-life and simulated series. Department of Statistical Sciences, University of Padua, Working Paper Series, n. 14. (http://www.stat.unipd.it/ricerca/fulltext?wp=422)
- Di Fonzo, T., Marini, M. (2011a), “Simultaneous and Two-step Reconciliation of Systems of Time Series: Methodological and Practical Issues”, Journal of the Royal Statistical Society C, 60, 143-164.
- Di Fonzo, T., Marini, M. (2011b), A Newton’s Method for Benchmarking Time Series according to a Growth Rates Preservation Principle. Department of Statistical Sciences, University of Padua, Working Paper Series, n. 7. (http://www.stat.unipd.it/ricerca/fulltext?wp=432)
- Dolan, E.D., Morè, J.J. (2002), “Benchmarking optimization software with performance profiles”, Mathematical Programming, 91, 201-213.
- Eurostat (2007), “Demetra. Advanced Seasonal Adjustment Interface”. (http://circa.europa.eu/irc/dsis/eurosam/info/data/demetra.htm)
- Fagan, J. (1995), Benchmarking. Unpublished note (mimeo).
- Fletcher, R., Reeves, J. (1964), “Function minimization by conjugate gradients”, The Computer Journal, 7, 149-154.
- Fletcher, R. (1987), Practical methods of optimization vol. 1: Unconstrained optimization. Wiley, New York.
- Fortier, S., Quenneville, B. (2009), Reconciliation and Balancing of Accounts and Time Series. From Concepts to a SAS procedure. In JSM 2009 Proceedings, Business and Economic Statistics Section. Alexandria, VA: American Statistical Association, pp. 130-144.
- Gilbert, J.C., Nocedal, J.(1992), “Global convergence properties of conjugate gradient methods for optimization”, SIAM Journal on Optimization, 2, 21-42.
- Griva, I., Nash, S.G., Sofer, A. (2009), Linear and Nonlinear Optimization. SIAM, Philadelphia.
- Hager, W.H., Zhang, H. (2006), “A survey of nonlinear conjugate gradient methods”, Pacific Journal of Optimization, 2, 35-58.
- Harvill Hood, C.C. (2005), An empirical comparison of methods for benchmarking seasonally adjusted series to annual totals. Paper presented at the Workshop on Frontiers in Benchmarking Techniques and their Applications to Official Statistics, Luxembourg, 7-8 April 2005.
- Hestenes, M.R., Stiefel, E.L. (1952), “Methods of conjugate gradients for solving linear systems”, Journal of Research of the National Bureau of Standards, 49, 409-436.
- Mittelmann, H.D., Pruessner, A. (2006), “A server for automated performance analysis of benchmarking data”, Optimization Methods and Software, 21, 105-120.
- Monsour, N.J., Trager, M.L. (1979), Revision and benchmarking of business time series. American Statistical Association, Proceedings of the Business and Economic Statistics Section, 333-337.
- Nocedal, J., Wright, S. (2006), Numerical Optimization. 2nd edition. Springer, New York.
- Polak E., Ribière G. (1969), “Note sur la convergence de directions conjugèes”, Revue Francaise d'Informatique et Recherche Opèrationelle, 16, 35-43.
- Polyak T. (1969), “The conjugate gradient method in extremal problems”. U.S.S.R. Computational Mathematics and Mathematical Physics, 9, 94-112.
- Powell, M.J.D. (1977), “Restart procedures for the conjugate gradient method”, Mathematical Programming, 12, 241-254.
- Powell, M.J.D. (1984), “Nonconvex minimization calculations and the conjugate gradient method”. In Lecture Notes in Mathematics, vol. 1066, Springer, Berlin, 122-141.
- Quenneville, B., Huot, G., Cholette, P., Chiu, K., Di Fonzo, T. (2003), Adjustment of seasonally adjusted series to annual totals. In Proceedings of Statistics Canada Symposium 2003 - Challenges in Survey Taking for the Next Decade, 2-12.
- Schmidt, M. (2006), The minFunc Toolbox for Matlab. (http://www.cs.ubc.ca/~schmidtm)
- The MathWorks (2009), Optimization Toolbox 4 User's Guide. Natick, MA.
- Titova, N, Findley, D., Monsell, B.C. (2010), Comparing the Causey-Trager method to the multiplicative Cholette-Dagum regression-based method of benchmarking sub-annual data to annual benchmarks. In JSM 2010 Proceedings, Business and Economic Statistics Section, 3007-3021.
- Trager, M.L. (1982), Derivation of Solution to the Benchmarking Problem: Relative Revision. Unpublished research notes, U.S. Census Bureau, Washington D.C.. Available as an appendix in Bozik and Otto (1988).
- U.S. Census Bureau (2009), X-12-ARIMA Reference manual, Version 0.3. U.S. Census Bureau, U.S. Department of Commerce, Washington D.C.. (http://www.census.gov/srd/www/x12a/)

### Appendix A. Feasible direction according to Causey and Trager (1981)
- Problem statement:
  - Given an (n x 1) (unconstrained) direction vector d, find an (n x 1) vector v solution to:
    - min. .1  and   = TT st v vdvvAv 0 . (29)
  - Desired properties: unitary norm (1 T vv), and feasibility (Av 0).
- Lagrangean function (equation (30)):
  -  * 1() 2 TTT L    vdvvAv
- First order conditions summarized (equation (31)):
  - 0 = L/v, 0 = L/dvT, 0 = L/vvT, 0 = L/Av
  - Presented compactly as:
    - 0 v
    - dvA 0
    - vv Av0
- Solution steps and intermediate results:
  - Pre-multiplying the first equation in (31) by A and solving for  yields (32):
    - 1 () T   AAAd
  - Substitution gives (33):
    - *1 () TT n       vIAAAAd
  - Noting N = 1 () TT n     A A A I is idempotent (T NNN) and T vv 1, squaring both sides of (33) yields:
    -  1 *2* 2 () TTT    vv dNddNd
  - Denoting 1  VN and substituting in (33) produces the final stated result:
    - vVd

### Appendix B. Performance Profiles (figures and captions)
- Figures present performance profiles for benchmarking algorithms on two series sets (EUQSA and MRTS) with two accuracy thresholds:
  - Acceptable solutions: 1 10  
  - Very accurate solutions: 3 10  
- Plotted performance metrics (for each figure):
  - Runtime (log2()) vs Percent of benchmarked series
  - Iterations (log2()) vs Percent of benchmarked series
  - Function evaluations (log2()) vs Percent of benchmarked series
- Algorithms compared (legend entries shown in figures):
  - Newton
  - Steepest Descent
  - Conjugate Gradient
  - Quasi-Newton BFGS
  - BMK1

*Source: _wp11179 - REFERENCES*

---


_Source: https://www.imf.org/-/media/websites/imf/imported-full-text-pdf/external/pubs/ft/wp/2011/_wp11179.pdf_