Thanh Viet Nguyen
email: thanhmpa@gmail.com

Michel Simioni

Dao Le Van

Assessment of TFP change at provincial level in Vietnam: New evidence using Färe-Primont productivity index

Keywords: Total factor productivity, Technical change, Technical efficiency, Mix and scale efficiencies, Färe-Primont index, Vietnam JEL Classification: D24, B41, B21, F63, O47, O53

Vietnam has become a lower middle-income country in less than 30 years, and is now facing the middle-income trap risk. Knowledge of changes in total factor productivity (TFP) is an essential element in assessing this risk. An in-depth analysis of the evolution of TFP and its determinants in Vietnam is presented in this paper. TFP evaluation uses a recently

Postprint of the article: Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis andPolicy, 64, 329-345. DOI : 10.1016/j.eap.2019.09.007 proposed multiplicative-complete economically ideal index, namely the Färe-Primont index, to evaluate TFP and to decompose it into its different components: technical change, pure technical, mix and scale efficiencies. TFP is computed at the provincial level over the 2010-2017 period. The results shows that estimated provincial TFP values are, on average, small whatever the considered year, but they have increased with an annual compound growth rate of 3.46%. Technical progress as measured by TFP* appears to be the main driver of TFP growth over the period, with an annual compound growth rate of 3.34%.

The expansion of the production set under constant returns-to-scale, from which TFP* is measured, is guided by movements of Ho Chi Minh city. Accordingly, on average, overall productive efficiency stagnated, with an annual compound growth rate of 0.12%. Technical efficiency has also stagnated over the period with its annual compound growth rate -0.62%.

The results imply that there has been an increasing gap between provinces in terms of the resource allocation efficiency. This evolution may have negative consequences on sustainable economic development and lead the country into the risk of middle income trap in the future.

Introduction

Vietnam's economic growth has been spectacular in the last three decades, changing the country from one of the world's poorest with Gross National Income (GNI) per capita around 527 constant 2010 US dollars in 1994 to a lower-middle income country with GNI per capita at 1741 constant 2010 US dollars in 2017 [START_REF] Fantom | The World Bank's classification of countries by income[END_REF]. But the foundations of Vietnamese growth are still fragile. [START_REF] Le | Winners and losers in Vietnam equitisation programs[END_REF] mentioned that since the introduction of "Doi Moi"policy in an attempt to move Vietnam towards a market economy, the transformation process has been slow and incomplete due to the remaining heavy influence of policies and institutions from the central planning days. Recent papers have pointed out that Vietnam could fail to transition to a high-income economy due to rising costs and declining competitiveness [START_REF] Pincus | Why doesn't Vietnam grow faster?: State fragmentation and the limits of vent for surplus growth[END_REF][START_REF] Herr | The integration of Vietnam in the global economy and its effects for Vietnamese economic development[END_REF][START_REF] Ohno | The quality of industrial policy and middle income traps: Comparing Vietnam with other countries[END_REF]. These papers discuss the risk of "middle income trap"faced by Vietnam. [START_REF] Ohno | The quality of industrial policy and middle income traps: Comparing Vietnam with other countries[END_REF] defines a middle income trap as "a situation where an economy is unable to create new value beyond what is delivered by given advantages". Given advantages include natural, demographic and geographical factors as well as external factors such as trade, aid and foreign investment. Development in the true sense occurs when value -added (GDP) is created and constantly augmented by domestic citizens and enterprises. Growth in Vietnam has been largely dominated by foreign-owned firms, and economic liberalization has been successful in making Vietnam regionally and globally integrated [START_REF] Ohno | The quality of industrial policy and middle income traps: Comparing Vietnam with other countries[END_REF]. But, such growth engine could sputter and lose power one day. As emphasized by [START_REF] Herr | The integration of Vietnam in the global economy and its effects for Vietnamese economic development[END_REF], "if the country does not manage to increase productivity permanently and innovative power, and at the same time create sufficient aggregate demand to keep the economy growing, a middle income trap becomes likely."

Recent analytical and empirical literature on middle-income traps has been surveyed by Agenor

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007 (2017). This paper argued that middle-income countries may end up being caught between lowwage poor countries, dominant in mature industries, and innovative rich countries, dominant in technology-intensive industries. [START_REF] Eichengreen | When fast growing economies slow down: international evidence and implications for China[END_REF] proposed an analysis of country characteristics and circumstances on which the timing of growth slowdown in fast-growing middle income countries depends. They found that around 85% of the growth slowdown is explained by the decrease in the TFP growth. More evidence can be found in [START_REF] Bulman | Transitioning from low-income growth to highincome growth: is there a middle-income trap?[END_REF] which showed that countries that managed to successfully overcome the middle-income range had relatively high TFP growth. [START_REF] Tho | The Middle-Income Trap: Issues for Members of the Association of Southeast Asian Nations[END_REF] claimed that middle income countries have to complete the "transition from input-driven to TFP-driven growth." The success stories of East Asia was supported by strong TFP growth, especially in China and Taiwan Province of China, where TFP contributed for more than half of all GDP per capita growth [START_REF] Aiyar | Growth Slowdowns and the Middle-Income Trap[END_REF].

A better knowledge of the evolution of productivity and its determinants in a country is a prerequisite for assessing the middle-income trap risk faced by this country. Various works give figures for Vietnam. According to estimates made by Vietnam National Productivity Institute, TFP accounted for about 48.5% in 2015 to Vietnam's economic growth and for over 30% in 2011-2015(VNPI, 2015)). [START_REF] Barker | Vietnam: The next Asian tiger? Economics discussion papers no[END_REF] showed also that labor productivity improvements accounted for 83% of the average growth by 5% of GDP per capita over the 1986-2014 period.

Vietnam's labor productivity has tripled from 2000 to 2017, and the gap with other comparable countries has narrowed [START_REF] Vnpi | Vietnam Productivity Report 2017[END_REF]. However, it should be noticed that Vietnam has a high proportion of agricultural workers (one half of total employment in 2013), and so, productivity in this country is still low. Indeed, productivity in the agricultural sector is generally lower than that in the industrial or service sectors. For instance, Singapore's labor productivity was 21 times higher than that in Vietnam In 1990, but only 12 times in 2016 (VNPI, 2015[START_REF] Vnpi | Vietnam Productivity Report 2017[END_REF].

This paper aims to contribute to the literature on middle-income trap risk in Vietnam, by providing a deeper evaluation of total factor productivity and its evolution using data on the 63

Vietnamese provinces over the 2010-2017 period. This contribution is threefold. First, this paper differs from existing literature which focuses primarily on labor productivity, by evaluating total

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019) factor productivity in a multiple outputs -multiple inputs framework. Technology is specified by disaggregating total provincial production in three components: agriculture, manufacturing, and services and considering three inputs: labor, capital and land. Second, this paper makes use of Färe-Primont index in order to measure total factor productivity. This index, which was introduced by O'Donnell (2014), belongs to the family of "multiplicative-complete economically ideal indices." These indices comply to all economically relevant axioms and tests defined by index number theory. Especially, the Färe-Primont index fulfills the identity axiom and the transitivity test, while the most commonly used productivity index, i.e. Malmquist index, fails to satisfy these properties (O'Donnell, 2012a). Total evolution of productivity over the studied period can then be decomposed in its evolution over smaller periods in a consistent way using Färe-Primont index.

Third, total factor productivity can be easily decomposed in its main drivers, i.e., technical change, pure technical efficiency change, mix efficiency change and scale efficiency change, using various Data Envelopment Analysis (DEA) linear programs. Special attention can then be devoted to the evolution of these productivity drivers not only over the entire period, but also over sub-periods.

Moreover, it is possible to characterize whether the Vietnamese provinces have evolved differently and to see if there are gaps between them, drawing policy implications at their disaggregated level instead than at only the national level.

The article is organized as follows. Section 2 provides an overview of existing literature on total factor productivity in Vietnam. Section 3 presents Färe-Primont productivity index and its decomposition into a measure of technical change and various measures of efficiency change including pure technical efficiency change, mix efficiency change and scale efficiency change.

Section 4 gives a description of the data. Section 5 is devoted to results presentation. Section 6 draws some policy implications for sustainable growth in Vietnam.

Literature Review

Research has been devoted to the impact of factor productivity on growth in Vietnam. These assessments are essentially made at the macroeconomic level. For instance, [START_REF] Park | Total factor productivity growth for 12 asian economies: The past and the future[END_REF] studies

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019) We conclude this overview of macroeconomic works on the impact of TFP on Vietnam's economic growth, mentioning the very recent work presented in [START_REF] Barker | Vietnam: The next Asian tiger? Economics discussion papers no[END_REF] [START_REF] Good | Index number and factor demand approaches to the estimation of productivity[END_REF], this paper

shows that firm productivity level increased after trade liberalization that occurred in 2007 when

Vietnam joined the World Trade Organization. Moreover, resource reallocation between firms was facilitated after the liberalization. [START_REF] Nguyen | Business reforms and total factor productivity in Vietnamese manufacturing[END_REF] shows also that Vietnamese firm-productivity increased over the 2000-2010 period, using also GSO data and applying a semiparametric method proposed by [START_REF] Wooldridge | On estimating firm-level production functions using proxy variables to control for unobservables[END_REF] and [START_REF] Petrin | Measuring aggregate productivity growth using plant-level data[END_REF] to measure firm-level TFP.

However, this evolution was contrasted according to sectors and regions. [START_REF] Huang | A new approach to estimating the metafrontier production function based on a stochastic frontier framework[END_REF]. This paper shows that, on average, SMEs can increase their current outputs by eight percent using the same quantity of inputs. Firms operating in major cities such as Hanoi and Ho Chi Minh City are found to be more efficient and possess better technology. Results indicate also that most SMEs in Vietnam use relatively lowlevel technologies, evidenced by the higher return from labour and raw materials than that from capital.

Our paper is halfway between these two literatures. Indeed, it is based on disaggregated data at the level of the provinces of Vietnam. But, unlike the macroeconomic or microeconomic works cited above, which are based on the assumption of single-product technology, it proposes a disaggregation of output into three components: agriculture, manufacturing and services. Computation of total factor productivity is not based on either a purely accounting approach or parametric assumptions about technology such as Cobb-Douglas (see the discussion of this assumption in [START_REF] Thai | Aggregate factor productivity through the ghosh approach: The paradox of the vietnamese economy[END_REF]. Our paper makes use of recent advances on TFP computation using multiplicative-complete economically ideal indices. These indices have good properties, including that of transitivity, and allow for a consistent assessment of the evolution of provincial TFP year by year without strong assumptions such as in previous papers. 2017) give the complete list of published empirical applications which includes [START_REF] Baležentis | The sources of the total factor productivity growth in Lithuanian family farms: a Färe-Primont index approach[END_REF], [START_REF] Islam | Broadacre farm productivity and profitability in south-western Australia[END_REF], [START_REF] Khan | Nonparametric estimates of productivity and efficiency change in Australian Broadacre Agriculture[END_REF][START_REF] O'donnell | Econometric estimation of distance functions and associated measures of productivity and efficiency change[END_REF], [START_REF] Rahman | Six decades of total factor productivity change and sources of growth in Bangladesh agriculture (1948-2008)[END_REF], and [START_REF] Tozer | Decomposing productivity and efficiency among Western Australian grain producers[END_REF] for agriculture; [START_REF] Widodo | Agglomeration economies and productivity growth in manufacturing industry: empirical evidence from Indonesia[END_REF] for manufacturing industry; [START_REF] Laurenceson | New estimates and a decomposition of provincial productivity change in China[END_REF] for provinces of China; [START_REF] Nguyen | Productivity and efficiency of Vietnamese banking system: new evidence using Färe-Primont index analysis[END_REF] for Vietnamese banks; and, [START_REF] Färe | Productivity change and fleet restructuring after transition to individual transferable quota management[END_REF] for fishery activities. See also [START_REF] Kar | Changes in total factor productivity and efficiency of microfinance institutions in the developing world: A non-parametric approach[END_REF] on microfinance institutions.

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007

posed by O'Donnell (2012a) and O'Donnell (2012b). These papers introduced a general class of multiplicatively-complete TFP indexes. The TFP index is defined as the ratio of an aggregate output to an aggregate input, and the change in TFP can then be expressed as the ratio of an output quantity index to an input quantity index, i.e. a measure of output growth divided by a measure of input growth. This means that, for province n in period t, TFP is given by

T F P nt = Y nt X nt (1)
where Y nt = Y (y nt ) and X nt = X(x nt ) represent the aggregate output and input, respectively, with y nt and x nt being the output and input vectors, respectively, and Y (.) and X(.) the aggregator functions.

Different aggregator functions give rise to different TFP indexes. A detailed list of usual aggregator functions, among them we find the usual Paasche and Laspeyres indexes, is given in O'Donnell (2012a). Among all the corresponding TFP indexes, we choose to compute the Färe-Primont index defined by O'Donnell ( 2014). Aggregator functions Y (.) and X(.) for this index are defined as

Y (y) = D O (x 0 , y, t 0 ) and X(x) = D I (x, y 0 , t 0 ) (2)
where

D O (x, y, t) = min{p > 0 :
x can produce y/p in period t} and D I (x, y, t) = max{p > 0 : x/p can produce y in period t} are, respectively, the Shephard output and input distance functions representing the technology available at period t, and x 0 and y 0 are, respectively, reference values of input and output for a representative time period t 0 [START_REF] Shephard | Theory of Cost and Production Function[END_REF].

In practice, the Färe-Primont index should be evaluated by choosing reference values that

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007

are relevant to the observations that are compared. For instance, if comparisons are to be made between all T observations in the data set, then possible choices for the reference values are the average quantities of outputs and inputs for each province computed over the observed period, i.e. x 0 = {x 0i } N i=1 and y 0 = {y 0i } N i=1 , with x 0i = T t=1 x it /T and y 0i = T t=1 y it /T . The representative period corresponds then to an hypothetical sample of provinces producing their sample average output quantities using their sample average input quantities. Then, DEA methodology can be used to compute the distance functions involved in the definition of Färe-Primont index, i.e. Eq. ( 2).

The Färe-Primont index can be shown as multiplicative-complete economically ideal in the sense that it satisfies all economically relevant axioms and tests from the index number theory: identity, transitivity, circularity, homogeneity, proportionality, time-space reversal and weak monotonicity axioms (see O'Donnell, 2012a). Moreover, unlike indexes such as Paasche and Laspeyres whose computation requires not only input and output quantities but also input and output prices, the computation of Färe-Primont index only requires observation of the quantities, not of the prices, which will be the case in our application.

Decompositions of TFP change O'Donnell (2012a) and O'Donnell (2012b) showed that all multiplicatively complete indexes can be decomposed into a measure of technical change and various measures of efficiency change. They first showed that the overall productive efficiency of a province, or TFPE, can be measured as the ratio of observed TFP of the province to the maximum TFP that is possible using the technology available in the considered period. The overall productive efficiency of province n in period t is thus

T F P E n,t = T F P nt T F P * nt = Y nt /X nt Y * t /X * t (3)
where T F P * t denotes the maximum achievable TFP using period-t technology, with X * t and Y * t denoting, respectively, the aggregate input and the aggregate output at this TFP-maximizing point.

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van Consider Fig. 1 where we report all potential combinations of aggregate output and input. Let point A represent this combination for a given province. Its TFP is measured by the slope of the line OA (the point O denoting the origin, with aggregate input and output quantities equal to zero). Let maximum achievable TFP* in the same period defined by the slope of the line OE.

Therefore, the overall productive efficiency of province A will be measured as the ratio of the slope of OA to the slope of OE. 3) can be decomposed into several ways using various efficiency measures. For instance, they define an output-oriented decomposition of the overall productive efficiency province n in period t as

T F P E n,t = OT E n,t × OM E n,t × ROSE n,t (4)
where OT E nt , OM E nt and ROSE nt denote measures of output-oriented pure technical efficiency, mix efficiency and residual scale efficiency, respectively.

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van The OTE measure is the well-known Farrell measure of technical efficiency [START_REF] Farrell | The measurement of productive efficiency[END_REF].

It measures pure technical efficiency as it compares the aggregate output of the province to the maximum quantity of aggregate output it could have produced using the same amount of aggregate input, keeping fixed the proportion of each output in the mix of outputs. Put differently, let the curve passing through point B represent the output mix-invariant production frontier in Fig. 1. Then OTE measures the increase in TFP that occurs when the province moves from point A to point B on the mix-invariant frontier. Put differently, OTE = slope of OA/ slope of OB.

The OME is a measure of the increase of TFP that can be gained now by holding inputs fixed and relaxing restrictions on output mix. This gain is measured by the ratio of the slope of OB to the slope of OC in Fig. 1, where point C belongs to the unrestricted or true production frontier, i.e. the boundary of the production possibilities set when all mix restrictions are relaxed.

Any increase in technical and mix efficiencies implies a rise in province TFP. When a province moves from point A to point C in Fig. 1, it becomes technically efficient and mix efficient. The provice increases the amounts of outputs it produces from fixed inputs, not only increasing initially these quantities while keeping the proportions between them fixed, but also by changing in a second time the proportions between the outputs. But province TFP is not yet maximized.

Province TFP will only be maximized by moving to the point D in Fig. 1. This point belongs to the straight line through the origin O, which is tangential to the true production frontier. This point thus defines the maximum attainable productivity given the technology at the considered period, or TFP*, or, put differently, the true production frontier with constant returns-to-scale.

The difference between the TFP at points C and D is defined as the residual output-oriented scale efficiency measure, or ROSE. In other words, residual output-oriented scale efficiency is a measure of the difference between TFP at a technically and mix efficient point and TFP at the point of maximum attainable productivity. Mathematically, ROSE is the ratio of the slope of OC to the slope of OD, or, similarly, to the slope of OE,.

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van To sum up, it can be easily checked that

T F P E n,t = Slope of OA Slope of OE = Slope of OA Slope of OB × Slope of OB Slope of OC × Slope of OC Slope of OE (5)
and that, by construction, each component in this decomposition is smaller or equal to 1.

The decomposition in Eq.( 4) focuses on the part of the efficiency of the province coming from a misallocation in the mix of outputs, and scale efficiency appears then as a residual. An alternative decomposition is also possible, namely

T F P E n,t = OT E n,t × OSE n,t × RM E n,t (6)
where OT E nt , OSE nt and RM E nt denote now measures of output-oriented pure technical efficiency, scale efficiency and residual mix efficiency, respectively.

The OTE measure still has the same interpretation in terms of pure technical efficiency. But now the decomposition focuses on scale efficiency, mix efficiency appearing only as a residual. Indeed, OSE now measures the gain in TFP a province can achieve by moving from the mix-invariant production frontier to the corresponding constant returns-to-scale mix-invariant production frontier. Output quantities, and thus aggregate output, will increase in order to reach the straight line, which is tangential to the mix-invariant production frontier (line OF in Fig. 2), i.e. achieving constant returns-to-scale but holding constant the output mix. Mathematically, OSE is equal to the ratio of the slope of OB to the slope of OG. Finally, the difference between the points G and D in Fig. 2 measures the residual part due to misallocation in the output mix. Indeed, we compare aggregate outputs belonging to two production frontiers with constant returns-to-scale but differing by the assumption on the invariance or not of the output mix. We denote by RME3 this measure whose value is equal to the ratio of the slope of OG to the slope of OD, or, similarly,

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van The last two terms in the previous two decompositions give the same value, which we denote by OSME for output-oriented mix and scale efficiency, i.e.

OSM E n,t = OM E n,t × ROSE n,t = OSE n,t × RM E n,t (8)
Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van TFP change for a given province n between two periods t 1 and t 2 can then be defined as

T F P n,t 1 ,t 2 = T F P n,t 2 T F P n,t 1 (10)
and, using the previous decompositions in Eqs. ( 4) and ( 6), decomposed as either

T F P n,t 1 ,t 2 = ( T F P * t 2 T F P * t 1 ) × ( OT E n,t 2 OT E n,t 1 ) × ( OM E n,t 2 OM E n,t 1 ) × ( ROSE n,t 2 ROSE n,t 1 ) (11)
or

T F P n,t 1 ,t 2 = ( T F P * t 2 T F P * t 1 ) × ( OT E n,t 2 OT E n,t 1 ) × ( OSE n,t 2 OSE n,t 1 ) × ( RM E n,t 2 RM E n,t 1 ) (12)
The first and second ratios in Eqs. ( 11) and ( 12), denoted by dTECH and dOTE, respectively, measure technological change and pure output-oriented technical efficiency change between the two periods. The last two ratios in Eq. ( 11), denoted by dOME and dROSE, respectively, measure changes in mix efficiency and in residual scale efficiency, respectively whereas the last two ratios in Eq. ( 12), which we denote by dOSE and dRME, respectively, measure changes in scale efficiency and residual mix efficiency, respectively. A value of dTECH larger than 1 indicates technical progress and smaller than 1 technical regress. The other ratios are efficiency changes with values larger than 1 indicating more efficiency and smaller than 1 indicating less efficiency relative to reference technologies in periods t 1 and t 2 .

The decompositions given in Eq. ( 11) and ( 12) will allow us to identify the main sources of productivity changes for each Vietnamese province.

In this research context, inputs must include the three main resources for growth and development at the provincial level, i.e. labor, capital and land. Outputs include total values of production in the following three sectors: agriculture, industry, and services.

Inputs Labor is measured by official number of workers aged over 15 in a province from General Statistical Office of Vietnam (GSO, 2015[START_REF] Gso | Yearbook statistics 2016[END_REF](GSO, , 2017(GSO, , 2018)). This measure is known having some limitations. First, it does not include half-time workers that may be present in agriculture, and self-employed workers. Second, some activities use workers under 15 years of age, which are not recorded too. However, despite these limitations, the official number of workers aged over 15 is considered as the best figure capturing labor force in Vietnam up to now.

General Statistical Office of Vietnam does not provide a measure of capital stock. Only data on total investment at provincial level are provided. It is well known that total investment is only a small amount of capital stock to be analyzed. Nevertheless, total investment can be used to recover capital stock using perpetual inventory method (OECD, 2009). Capital at time t is thus defined as

K t = (1 -δ) × K t-1 + I t , t = 1, . . . , T (13)
where K t denotes capital stock in year t, I t , total investment in year t, and δ, depreciation rate.

In Eq. ( 13), the total investment series is known, but not the capital series. This latter can be initialized using K 0 = I 0 /(δ + θ) where I 0 denotes total investment of the initial year, and θ is the growth rate of gross output over the period, computed as θ = (GDP T /GDP 0 ) 1/T where T is the last year of observation. Depreciation rate δ is computed as the average of depreciation rates over the studied period.

Land is considered as an input. Data on agricultural land, non-agricultural land and other lands provided by General Statistical Office of Vietnam (GSO, 2015[START_REF] Gso | Yearbook statistics 2016[END_REF](GSO, , 2017(GSO, , 2018) ) Outputs Most papers dealing with productivity measurement at the provincial level use provincial Gross Domestic Product (GDP) as output. Hereafter, GDP is divided into the total value of products in agriculture, the total value of products in industry and total value of products in services. All data are converted into 2010 Vietnamese Dong for ease of comparison and evaluation.

Table 1 lists some descriptive statistics of the input and output data.

Results Analysis

Total factor productivity and technical change The results from calculating TFP and decomposing it in its main components, as explained in section 3, are shown in The fact that, although growing, TFP levels remain low, while we observe technical progress at the same time, can be explained by considering the evolution of the set of attainable productions combinations, or production set, between 2010 and 2017. The top panel in Figure 3 reports the estimated production frontiers using aggregated input-output combinations in 2010 and 2017, and the location of each province these two years. 5 The production set has expanded over the period and this expansion appears to be largely determined by TFP evolution of Ho Chi Minh City province as indicated by its location in 2010 (HCM2010) and 2017 (HCM2017). But, the comparison of other provinces locations does not show any major evolutions of them in the production set compared to Ho Chi Minh City, except for some of them, including Hanoi.

Bottom panel in Figure 3 shows more clearly some of these changes and their location in relation to the production frontier. The observed expansion of the production frontier is guided by

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019) The availability of TFP estimates for each province over the eight years period also allows to see how its spatial distribution has evolved. one. The high and stable TFP growth of these regions have been documented in several studies (see [START_REF] Nguyen | Business reforms and total factor productivity in Vietnamese manufacturing[END_REF][START_REF] Le | Technical efficiency of small and medium manufacturing firms in vietnam: A stochastic meta-frontier analysis[END_REF]. Central Highlands region had known the fastest TFP growth rate (with the greatest color shift). The Central Highlands region is also increasingly contributing to the sustainable growth of Vietnam in the period from 2010 to now.

Technical efficiency and scale-mix efficiency TFPE can be decomposed into OTE and OSME components as shown in section 3. An estimated value of OTE equal to one indicates that the corresponding province is located on the boundary of the output mix-invariant production set, and thus is technically efficient given the mix of outputs this province produce. An estimated value below one means that the province is located under the output mix-invariant production frontier

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007 and hence is technically inefficient. That is, given the input quantities used, output quantities can then be proportionally increased, or keeping fixed the mix of these outputs. Similarly, an estimated value of OSME equal to one indicates that, given the input quantities used, the province produce the maximum attainable aggregate output quantity. This province has chosen the optimal mix of outputs and has exhausted any scale-economies. An estimated value below one means that the province has not fulfilled this two objectives.

Average value of OTE and OSME over the 2010-2017 period are reported in Table 3. The estimated average values of the two indicators are very far from one. On average, Vietnamese provinces exhibited technical, and scale and mix inefficiencies. In addition, both indicators have experienced contrasting trends, which explains the stagnation of their product, namely TFPE.

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007 At the same time, scale and mix efficiency has grown from 57.53% to 60.36%. OSME indicator combining two potential sources of inefficiency is difficult to interpret as it stands. We will return to the interpretation of the evolution of this indicator during its decomposition below.

The results given in the Table 3 are average results and it may be interesting to look in more detail at the evolution of the distributions of OTE and OSME as given in the Yen Bai, Tuyen Quang, Bac Kan, Lang Son, Phu Tho and some provinces scattered across the country. These are provinces with a relatively small size of GDP and a greater potential for more efficient resource reallocation. In contrast, Ho Chi Minh City and Hanoi have highest GDP values in the country and may encounter more difficulty when improving the resource allocation efficiency (input mix and output mix). Thus, OSME growth will make a significant contribution to TFP growth towards sustainable development, particularly in provinces with high GDP levels.

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019) which TFP* is measured, is guided by movements of Ho Chi Minh province where is located the Vietnam's economic capital city. Accordingly, on average, overall productive efficiency stagnated, with an annual compound growth rate of 0.12%.

• Stagnation of overall productive efficiency can be explained by the evolution of its components. There was a decrease in average technical efficiency over the period from 2010 to 2015 and it has increased in 2016 and 2017. Overall, technical efficiency stagnated over the period (its annual compound growth rate is only -0.62%). At the same time, average scale and mix efficiency increased but with also at very small compound growth rate of 0.83%. It is also noteworthy that significant gains can be made by decreasing both the two sources of inefficiency (23% and 40% for technical efficiency and scale and mix efficiency in 2017, respectively). When only focusing on either mix efficiency or scale efficiency, gains achieved by improving efficiency are around 11%.

• There is the existence of two groups of provinces when looking at OTE distributions. The first group, which increases in size during the period, is made up of relatively technically efficient provinces, i.e. with OTE score equal or very close to one. A second group, which decreases in size during the period, consists of provinces that are much less technically efficient, their score being between 0.4 and 0.8.

From our results, we can draw several policy implications. First, the government should pay more attention on R&D activities in order to improve technical progress for all provinces. In

2

  Färe-Primont productivity index For the purpose of this article, we use the recent developments in TFP index measurement and TFP index decomposition pro-According to Molinos-Senante et al. (2017), Färe-Primont index has been scarcely applied empirically. Molinos-Senante et al. (

Figure 1 :

 1 Figure 1: TFP definition and decomposition focusing on mix efficiency

Figure 2 :

 2 Figure 2: TFP definition and decomposition focusing on scale efficiency

  Figure 4 summarizes this evolution by considering three different years: 2010, 2013 and 2017. Maps show that Red River Delta and South-East regions have known the highest TFP in the country (darkest color). These regions are two bestdeveloped areas of the country in terms of natural conditions, capital, and location. Hanoi is located in Red River delta region while Ho Chi Minh City and Vung Tau are in the South-East

Figure 3 :

 3 Figure 3: The change in the production frontier over the period 2010-2017

Figure 4 :

 4 Figure 4: TFP change across Vietnam for 2010-2017

  Figure 5. The most striking result is the existence of two groups of provinces when looking at OTE distributions. The first group, which increases in size during the period, is made up of relatively technically efficient provinces, i.e. with OTE score equal or very close to one. A second group, which decreases in size during the period, consists of provinces that are much less technically efficient, their score being between 0.4 and 0.8. Meanwhile, the distribution of OSME is unimodal and its mode, Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007 around 0.4, has not changed during the period. But, this distribution becomes more flattened,with more and more provinces having an OSME value higher than the mode.

Figure 5 :

 5 Figure 5: Changes in TFPE, OTE, and OSME distributions over the period 2010-2017

Figure 6 :

 6 Figure 6: Changes in OSME, OME, and ROSE distributions over 2010-2017 period

Figure 7 :

 7 Figure 7: Changes in OSME, OSE, and RME distributions over 2010-2017 period

Figure 8 :

 8 Figure 8: Evolution of TFP and its components: OTE and OSME

Figure 9 :

 9 Figure 9: Evolution of OSME and its components: first decomposition

Figure 10 :

 10 Figure 10: Evolution of OSME and its components: second decomposition

  addition, it is important to have key economic regions to stimulate TFP growth in neighborhood thought technological learning. Second, technical efficiency need to be improved, especially in the Northwest, Northeast and South Central regions where technical efficiency has been low and decreased in size during the study period. Technical efficiency can be improved by (a) removing barriers to adopt new technologies, (b) providing education and training services to provincial decision makers about the existence and proper use of new technologies (e.g., agricultural extension programs), and (c) ensuring that output markets are competitive (e.g., by removing barriers to market entry). Finally, the gap of technical efficiency between two groups of provinces may have

  seven Asian countries, including Vietnam, and the impact of TFP on growth over the 1970-2007 period. Average TFP growth rate of these Asian countries is evaluated at 6.09% over this period, i.e. a higher rate than other regions in the world, using growth accounting model.

	Moreover, Park (2012) shows that TFP was only a minor contributor to growth over the 1970-
	2000 period and that the 2000-2007 period can be considered as transition toward productivity
	based growth. Using an econometric model of TFP growth, Park (2012) also forecasts that
	TFP will continue to increase in the Asian countries. In particular, TFP growth in Vietnam is
	forecasted to increase at a rate about 1.08% to 2.85% per year over 2010-2020 and about 1.09%
	to 2.82% over 2020-2030.
	More recently, VNPI (2015) provides and overview of labor productivity and growth evolutions
	in Vietnam over the 1990-2015 period. This study shows that labor productivity tripled from
	1990 to 2015, evolving from 2800 US Dollar (in terms of purchasing power parity) to 8400 US
	Dollar. Specifically, TFP growth contributed increasingly to GDP growth and TFP grew rapidly
	over the period 2011-2015. More precisely, Vietnam's TFP grew at an average annual rate of
	1.79% and contributed about 30% to GDP growth in this period. VEPR (2017) which study labor
	productivity and minimum wage contribution to economic growth for the 2009-2016 period, shares
	the same view on TFP growth and its contribution to Vietnamese economy growth. Moreover,
	VEPR (2017) shows that of excessive wage intervention policies have restricted growth potential
	in Vietnam.

. Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007 the growth of

  Northern, the Southern and the Central regions. Lastly,[START_REF] Le | Technical efficiency of small and medium manufacturing firms in vietnam: A stochastic meta-frontier analysis[END_REF] focuses on Small and Medium Enterprises (SME) in Vietnam. It aims at estimating technological gaps and identifying factors affecting variations in SMEs' technical efficiency using firm-level survey data in 2008 and stochastic meta-frontier framework of

	Most sectors have seen
	very limited growth, while the technology sector has the fastest growth rate. Moreover, firm
	productivity growth have been faster in the 2000-2005 period than in the 2005-2010 period.
	More sectors with positive and faster growth rate are observed in the 2000-2005 period in other
	areas rather than four key economic regions. 1 Slowdown of TFP growth is shown for several
	sectors in negative TFP growth rates in the 2005-2010 period, especially in other regions. The
	Southern key economic region, which is the biggest economic hub of Vietnam, performed at more
	stable TFP growth rate during the two periods. The youngest key economic region, i.e., Mekong
	Delta, and other areas were in deeper slowdown of TFP in the 2005-2010 period compared to the
	1 Key economic regions were assigned by the government since 1997 to take advantages of the local region's
	natural resources and comparative advantages as well as to support for satellite provinces. Four key economic
	regions in Vietnam are: (i) The Northern key economic region includes Ha Noi (capital), Hai Phong, Vinh Phuc,
	Bac Ninh, Hung Yen, Quang Ninh, and Hai Duong. (ii) The Central key economic region consists of Da Nang,
	Thua Thien Hue, Quang Nam, Quang Ngai, and Binh Dinh. (iii) The Mekong River Delta economic region covers
	the area of Can Tho, An Giang, Kien Giang, and Ca Mau. (iv) Provinces in the Southern economic region are
	Ho Chi Minh, Dong Nai, Ba Ria-Vung Tau,Binh Duong, Binh Phuoc, Tay Ninh, Long An, and Tien Giang.

Comment citer ce document : Nguyen, T. V., Simioni, M., Le

Van, D. (2019)

. Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007

  for fishing activities. Agricultural production land refers to land used in the agricultural production, including annual crop land and perennial crop land. Forestry land refers to land used in forestry production, including: productive forest, protective forest and specially used forest. Non-agricultural land includes special used land and homestead land. Specially used land is land being used for other purposes, not for agriculture, forestry and living. Homestead land is land used for housing and other construction works serving urban activities.

	are used
	to measure this input. Agricultural land includes agricultural production land, forestry land and

Comment citer ce document : Nguyen, T. V., Simioni, M., Le

Van, D. (2019)

. Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007 water surface land

Table 3

 3 Results in Table3show that average TFP has grown gradually from 0.3246 in 2010 to 0.4222 in 2017. This increase is mainly due to TFP*, which rose from 0.6887 in 2010 to 0.8977 in 2017, TFPE remaining relatively stable over the period. In other words, the observed growth in average total factor productivity over the 2010-2017 period stems solely from technical progress observed over the same period. Average overall factor productivity, meanwhile, remained stable at the same time.

. The first part of this table reports the geometric average values of provincial measures of total factor productivity, or TFP, maximum achievable total factor productivity, or TFP*, overall productive efficiency, or TFPE, pure output-oriented technical efficiency, or OTE, and output-oriented scalemix efficiency, or OSME. 4 Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007

Table 1 :

 1 Description of the outputs and inputs

		Inputs:			Outputs:
		Labor	Land	Capital	Agri. Industry Services
	Year = 2010
	Mean	776.54 525.49	140076 6583.73 15356.04 16343.66
	Std.dev.	77.57	46.39	26950.03	521.02	4115.26	4917.58
	Median	627.13 473.98	83893.26 5267.54	6266.33	8950.59
	Minimum	187.65	82.27	22987.85	975.00	889.13	1932.72
	Maximum	3705.62 1649.03 1232452.06 18768.30 209471.00 268991.00
	Year = 2011
	Mean	796.40 525.49	151105 6887.39 16819.27 17928.61
	Std.dev.	79.42	46.39	29124.87	554.72	4280.70	5458.90
	Median	647.60 473.98	91049.37 5535.77	6946.59	9986.89
	Minimum	195.89	82.27	27025.57	889.00	824.34	2084.98
	Maximum	3816.69 1649.03 1353072.72 20424.12	211565	302065
	Year = 2012
	Mean	812.48 525.50 162659.37 7217.52 18335.74 19585.11
	Std.dev.	81.68	46.39	31426.29	583.09	4663.71	5957.13
	Median	659.90 473.98	97527.35 5758.86	7659.91 10752.80
	Minimum	199.93	82.27	29677.37	862.00	801.55	2166.30
	Maximum	3943.18 1649.09 1466964.44 22617.10 229577.00 333275.00
	Year = 2013
	Mean	826.28 525.51 175700.81 7476.39 20123.42 21254.16
	Std.dev.	82.81	46.39	33881.66	608.00	4900.16	6530.79
	Median	670.20 473.98 103830.46 6082.93	8331.11 11470.70
	Minimum	203.19	82.27	31733.87	940.00	757.98	2334.27
	Maximum	3989.24 1649.00 1581740.04 23764.80 232546.00 368992.00
	Year = 2014
	Mean	834.55 525.44 190604.67 7781.14 21594.97 22895.31
	Std.dev.	83.80	46.37	36684.00	632.65	5088.96	6999.77
	Median	681.50 473.74 113066.33 6236.60	9512.30 12259.00
	Minimum	205.91	82.27	33199.04	879.00	755.06	2540.20
	Maximum	4059.16 1649.00 1701362.74 25263.80 238172.00 394734.00

Comment citer ce document : Nguyen, T. V., Simioni, M., Le

Van, D. (2019)

. Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007

Table 2 :

 2 Description of the outputs and inputs(cont'd)

		Inputs:			Outputs:
		Labor	Land	Capital	Agri. Industry Services
	Year = 2015
	Mean	841.05 525.83 206846.54 8069.10 23413.40 24778.70
	Std.dev.	84.81	46.35	39495.32	655.79	5262.60	7517.57
	Median	686.80 473.74 122906.55 6402.10	9852.45 13259.00
	Minimum	221.77	82.27	34076.40	848.00	827.88	2707.34
	Maximum	4129.54 1649.00 1837123.42 26317.70 236303.00 423243.00
	Year = 2016
	Mean	849.36 525.83 224021.38 8260.49 24706.96 26812.98
	Std.dev.	85.93	46.35	42717.80	663.94	5238.29	8126.82
	Median	690.81 473.74 133263.50 6413.00 10206.20 14384.30
	Minimum	226.98	82.27	35635.36	873.00	848.54	3050.40
	Maximum	4224	1649	1986521	25694	212136	457253
	Year = 2017
	Mean	860.05 525.81 242701.69 8512.26 27251.11 28977.07
	Std.dev.	87.32	46.29	46277.93	683.53	5521.52	8797.57
	Median	694.74 473.74 142091.14 6662.00 11609.00 15575.50
	Minimum	226.96	82.27	37271.52	900.00	893.14	3226.08
	Maximum	4320.61 1648.16 2150777.43 26439.10 203337.00 493996.00

  . Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007 the evolution of Ho Chi Minh City. This province achieve a high and stable TFP growth rate during 2010-2017 period. Hanoi also exhibits a high and stable TFP growth rate but slower than Ho Chi Minh City. Moreover, this growth is mainly driven by rapid growth in aggregate input. Growth of Hanoi seems unsustainable as aggregate output growth stayed proportional to aggregated input growth without a more efficient input use. Two other provinces, i.e. Thai Nguyen (TN) and Bac Ninh (BN), which are located in the north of Hanoi, have known similar TFP evolution as Hanoi although these evolutions are less pronounced. Indeed, these two provinces have known rapid capital growth through attracting large FDI. Finally, note the evolution of Vung Tau province (VN) whose TFP is closed the maximum attainable TFP level in each year and where input use is always efficient while not increasing. Vung Tau is the only petroleum base of Vietnam where crude oil and natural gas exploitation activities dominate province's economy and contribute to principal income to Vietnam's budget and export volume.

  Table 3 focuses on decomposing OSME into components in two ways. The first consists in decomposing OSME into OME and ROSE. This decomposition focuses on the impact of changing the output mix on efficiency. Average values of OME are quite stable and large, indicating that only small losses, between 9 and 12%, came from inefficient choice in the output mix, for fixed input use, over the period. The main source of inefficiency came then from residual scale efficiency, once output mix has been adjusted. The increase in ROSE average values brought about that of OSME over the period. citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007

	Comment

Table 4 :

 4 Table 4 lists the geometric average values of these measures. 6 An average value greater than unity indicates an improvement in the measure and a value less than unity indicates deterioration. Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007 Changes in total factor productivity and its component average growth rates between 3.27% and 4.95%. TFP* growth comes from conditions such as technological improvements, improved quality of resources, application of science in production management. Over the whole period TFP* growth potential reached 30.35%. The change in overall potential productivity, or TFPE, does not follow a clear pattern. TFPE has grown by only 0.92% over the period 2010-2017. Technical progress appears as the main driver of total productivity change. Pure technical efficiency, or OTE, decreased by 4.85% over the overall period. The change in OTE may have resulted in an increasing gap between provinces. This can have many negative consequences for socio-economic life, but is quite common in developing countries. Some regions with distinct advantages began to develop rapidly and thus created distances with other regions, as shown in China by[START_REF] Tan | Uneven growth of urban clusters in mega regions and its policy implications for new urbanization in China[END_REF] and[START_REF] Li | Uneven urban-region sprawl of China's megaregions and the spatial relevancy in a multi-scale approach[END_REF]. Meanwhile, the increase of OSME in the 2010-2017 period is 6.91% and contributed significantly to the growth of TFPE.The second part of Table4focuses on the decomposition of the change in OSME into changes in its components using either Eq. (11) or Eq. (12). In general, the change is small regardless of the focus of the analysis: change in output mix or change in scale. Indeed, OME decreased only

	The first part of table 4 shows that TFP changes throughout 2010-2017 were relatively stable,

Comment citer ce document : Nguyen, T. V., Simioni, M., Le

Van, D. (2019)

. Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007 with Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007

It can be easily verified that this measure is output-oriented as well as input-oriented, and hence the absence of O in its acronym

In this table and those that follow, we report geometric average values of individual Färe-Primont indexes and associated effciency measures. Indeed, by construction, Färe-Primont indexes are multiplicative and geometric mean has proved to be an adequate tool when measuring the central tendency of numbers whose values are meant to be multiplied together (see, for instance,[START_REF] Nguyen | Productivity and efficiency of Vietnamese banking system: new evidence using Färe-Primont index analysis[END_REF].

Estimates are recovered using classical DEA program under variable returns-to-scale assumption.

For each province, the changes are computed by taking the values of the different indicators for that same province in 2010 as references, using either Eq. (11) or Eq. (12).

Acknowledgements

This research is funded by Vietnam National University, Hanoi (VNU) under project number QG.17.35. Michel Simioni thanks INRA-CIRAD GloFoodS program for supporting his stay as visiting professor at Hanoi University for Natural Resources and Environment, Vietnam from October 2018 to February 2019.

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007 negative consequences on sustainable economic development for whole country and may lead the country into the risk of middle income trap in the future. Thus, it is necessary to improve the quality of education and training services as well as R&D activities for those provinces that are lagging behind so that they can catch up to those provinces on the production frontier.

Declaration of Competing Interest

All authors declare that they have no conflict of interest.

Comment citer ce document : Nguyen, T. V., Simioni, M., Le Van, D. (2019). Assessment of TFP change at provincial level in Vietnam: new evidence using Färe-Primont productivity index. Economic Analysis and Policy, 64, 329-345. , DOI : 10.1016/j.eap.2019.09.007


I Misztal

A Legarra

Invited review: efficient computation strategies in genomic selection

Keywords: genomic selection, single-step, genomic relationship matrix, inverse, REML Implications

The purpose of this study is review and evaluation of computing methods used in genomic selection for animal breeding. Commonly used models include SNP BLUP with extensions (BayesA, etc), genomic BLUP (GBLUP) and single-step GBLUP (ssGBLUP). These models are applied for genomewide association studies (GWAS), genomic prediction and parameter estimation. Solving methods include finite Cholesky decomposition possibly with a sparse implementation, and iterative Gauss-Seidel (GS) or preconditioned conjugate gradient (PCG), the last two methods possibly with iteration on data. Details are provided that can drastically decrease some computations. For SNP BLUP especially with sampling and large number of SNP, the only choice is GS with iteration on data and adjustment of residuals. If only solutions are required, PCG by iteration on data is a clear choice. A genomic relationship matrix (GRM) has limited dimensionality due to small effective population size, resulting in infinite number of generalized inverses of GRM for large genotyped populations. A specific inverse called APY requires only a small fraction of GRM, is sparse and can be computed and stored at a low cost for millions of animals. With APY inverse and PCG iteration, GBLUP and ssGBLUP can be applied to any population. Both tools can be applied to GWAS. When the system of equations is sparse but contains dense blocks, a recently developed package for sparse Cholesky decomposition and sparse inversion called YAMS has greatly improved performance over packages where such blocks were treated as sparse. With YAMS, GREML and possibly singlestep GREML can be applied to populations with >50 000 genotyped animals. From a computational perspective, genomic selection is becoming a mature methodology.

Introduction

The concept of genomic selection [START_REF] Meuwissen | Prediction of total genetic value using genome-wide dense marker maps[END_REF] generated great excitement in the animal breeding community. With genomic information from SNP panels, one can achieve accuracy from young animals almost as high as from a progeny test [START_REF] Schaeffer | Strategy for applying genome-wide selection in dairy cattle[END_REF][START_REF] Vanraden | Invited review: reliability of genomic predictions for North American Holstein bulls[END_REF] at a greatly reduced cost. Initially the genomic computations used versions of BLUP with SNP fitted as random effects, such as SNP BLUP, BayesB, etc., (see [START_REF] Gianola | Additive genetic variability and the Bayesian alphabet[END_REF] for a review). An alternative form of SNP BLUP is genomic BLUP (GBLUP), where the animal effect is fit with a GRM [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF]. For prediction, SNP BLUP and GBLUP are equivalent models [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Strandén | Allele coding in genomic evaluation[END_REF].

When only a small fraction of the population is genotyped, the information from ungenotyped animals can be summarized in pseudo-observations for genotyped animals. Alternatively, GBLUP can be extended to single-step GBLUP (ssGBLUP) (Aguilar et al., 2010; Christensen and Lund, 2010), where a numerator relationship matrix (NRM) for all individuals and GRM are combined and then applied to BLUP. Benefits of ssGBLUP include simplicity of application (another BLUP), avoidance of double counting, and accounting for pre-selection on Mendelian sampling [START_REF] Legarra | Single step, a general approach for genomic selection[END_REF]. ssGBLUP was † E-mail: Ignacy@uga.edu Animal, page 1 of 6 © The Animal Consortium 2016. This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution, and reproduction in any medium, provided the original work is properly cited. doi:10.1017/S1751731116002366 animal extended to be based partially or completely on SNP effects (Fernando et al., 2014; Liu et al., 2014).

Large interest in the genomic community exists in identifying an optimal set of SNP and their variances for increased accuracy of evaluation. Such an identification is straightforward with Bayesian SNP models (e.g. [START_REF] Gianola | Additive genetic variability and the Bayesian alphabet[END_REF]. As SNP BLUP and GBLUP are equivalent, SNP and their variances selected in Bayesian SNP models can be transformed to weighted GRM and subsequently BLUP. Methods exist to generate SNP weights directly in the GBLUP model [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Zhang | Best linear unbiased prediction of genomic breeding values using a trait-specific markerderived relationship matrix[END_REF][START_REF] Sun | A fast EM algorithm for BayesA-like prediction of genomic breeding values[END_REF] or ssGBLUP [START_REF] Wang | Genome-wide association mapping including phenotypes from relatives without genotypes[END_REF].

Initially the extent of genomic selection was limited by costs of genotyping. With these costs dropping dramatically, the number of genotyped animals has increased dramatically, with very high-density SNP chips available. In dairy cattle, over 1 million Holsteins had been genotyped in the United States as of March 2016 (Council on Dairy Cattle Breeding; https://www.cdcb.us/Genotype/cur_density.html). Subsequently, computing costs are becoming more of an issue. The purpose of this paper is to present and discuss computational methods in genomic selection.

Definition of terms and basic models

Let A denote a NRM based on pedigrees. Let Z be a matrix of gene content with z ij being the number of occurrences of the minor allele in SNP j of animal i. The values of 'raw' z ij are 0, 1 and 2, whereas that of centered z ij is -2p j , 1p j and 2p j , where p j is gene frequency of SNP j. The simplest SNP BLUP model is

y i = μ + X j z ij a j + e i
where y i is a (pseudo)phenotype of animal i, a $ N 0; Iσ 2 a À Á the vector of SNP effects, e $ N 0; Iσ 2 e À Á the vector of residuals, and σ 2 a and σ 2 e are SNP and residual variances, respectively. GBLUP can be defined as

y i = μ + u i + e i where u $ N 0; Gσ 2 u À
Á is a vector of animal effects, G the GRM and σ 2 u the additive variance. G can be derived by a transformation from SNP BLUP:

G = ZZ 0 σ 2 a σ 2 u
possibly with an alternative scaling factor [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF]:

G = ZZ 0 P j p j 1Àp j À Á Â Ã
If SNP BLUP involves heterogeneous variance for SNP: a $ Nð0; Dσ 2 a Þ, where D is a diagonal matrix of weights, the equivalent 'weighted' G is

G = ZDZ 0 P j p j 1Àp j À Á Â Ã
SNP BLUP and GBLUP are equivalent models [START_REF] Strandén | Allele coding in genomic evaluation[END_REF].

When only a fraction of animals is genotyped, the numerator and genomic relationships can be combined (Legarra et al., 2009; Christensen and Lund, 2010):

H = A + A 12 A À1 22 0 0 I ! I I ! ðGÀA 22 Þ I I ½ A À1 22 A 21 0 0 I !
where H is a combined matrix and indices in A refer to ungenotyped (1) and genotyped (2) animals. The inverse of matrix H is [START_REF] Aguilar | A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF][START_REF] Christensen | Genomic prediction when some animals are not genotyped[END_REF])

H À1 = A À1 + 0 0 0 G À1 ÀA À1 22 !
BLUP with matrix H is called ssGBLUP.

Computing operations in animals breeding

Common operations in animal breeding include solving of mixed model equations and variance component estimation either via REML or Bayesian methods with Gibbs sampling. Such operations in the animal breeding context are described in [START_REF] Mrode | Linear models for the prediction of animal breeding values, 2nd and 3rd edition[END_REF] and [START_REF] Misztal | Computational techniques in animal breeding[END_REF], and are typically based on Henderson's mixed model equations (MME). Let this linear system of equations be

B x = y
where B is left-hand-side (LHS), y the right-hand side, and x the vector of solutions. Solving mixed model equations can be accomplished by finite and iterative methods. In finite methods, exact solutions (except for numerical errors) are obtained in a finite number of steps. Popular methods for general matrices are based on the LU decomposition (B = LU, where L and U are lower and upper triangular, respectively) and for symmetric matrices the Cholesky decomposition (B = LL', where L is lower triangular). As often mixed model equations are not full rank, finite methods used for mixed models need to compute generalized solutions rather than cause numerical exceptions. For example, the Cholesky decomposition can be modified to

B = LDL′
where D is a diagonal matrix with zero elements for redundant equations and non-zeros otherwise. Although finite methods usually compute solutions with high precision, their cost is high, usually close to cubic. Finite methods can be used to calculate inverses, which are used in REML or for prediction error variance (PEV) calculations. In general, the cost of the inverse is a few times that of a solution by decomposition.

Misztal and Legarra

When matrices are sparse, with nearly all elements of LHS equal to 0, a 'sparse' Cholesky decomposition exists that avoids operations with zeros, leading to about quadratic rather than cubic costs (George and [START_REF] Liu | Computer solution of large sparse positive definite systems[END_REF]. Also, an algorithm exists for computing a 'sparse' inverse, where elements computed are those corresponding to non-zeroes in the original matrix [START_REF] Takahashi | Formation of a sparse bus impedance matrix and its application to short circuit study[END_REF]. Such an inverse is sufficient to calculate traces in REML [START_REF] Misztal | Sparse matrix inversion for restricted maximum likelihood estimation of variance components by expectationmaximization[END_REF] or PEV. However, equations involving SNP effects are usually dense. When the system of equations is sparse but contains dense blocks, a sparse matrix package can be modified for greatly improved performance in REML and other applications [START_REF] Masuda | Technical note: acceleration of sparse operations for average-information REML analyses with supernodal methods and sparse-storage refinements[END_REF]. This is quite new as most matrix packages are efficient for either dense or sparse matrices, but not for mixtures of both.

Iterative methods are such that progressively more accurate solutions are obtained every next round. Two such methods are popular in animal breeding. The first one is Gauss-Seidel iteration (GS), which is also a backbone of Gibbs sampling. GS is very simple, but requires access to equations by rows. The second one is preconditioned conjugate gradient (PCG), a method difficult to understand, but easy to implement and usually with superior convergence rate to GS [START_REF] Tsuruta | Use of the preconditioned conjugate gradient algorithm as a generic solver for mixed-model equations in animal breeding applications[END_REF].

In animal breeding application, the size of LHS can be much larger than the size of data to generate LHS. For example, one line of data in a 10 trait model with 10 effects per trait includes 20 numbers (one for each trait and effect) but generates 10 000 contributions to LHS. A special implementation of iteration algorithms is matrix-free or by iteration on data [START_REF] Schaeffer | Computing strategies for solving mixed model equations[END_REF][START_REF] Misztal | Indirect solution of mixed model equations[END_REF], where coefficients of LHS are regenerated from the data every round of iteration. Although the implementation of iteration on data by GS can be complicated as the coefficients need to be recreated row by row, such an implementation in PCG is trivial as PCG only requires a product of B by a vector (Bq, with q a given vector), with no individual elements necessary (except for diagonals).

Tricks and rules in mixed model and genomic computations

Product of numerator relationship matrix Computing the NRM A by a tabular method has a quadratic cost. Also, computing relationships for specific animals requires creating a matrix for all ancestors. For large populations, A would not fit into memory. When only a product A by a vector, say q, is needed, Aq can be computed inexpensively by two scans of the pedigree [START_REF] Colleau | An indirect approach to the extensive calculation of relationship coefficients[END_REF]. By selecting zeros in q, one can also compute a product of a section of A by a vector (e.g. A 22 q) (Misztal et al., 2009; Aguilar et al., 2011).

Sequential multiplication and preconditioned conjugate gradient

The LHS of the mixed model equations contains expressions like Z'Z. The PCG algorithm by iteration on data requires computing Z'Zq, where q is a given vector. [START_REF] Strandén | Solving large mixed linear models using preconditioned conjugate gradient iteration[END_REF] found that while the cost of (Z'Z)q is very high, the cost of Z'(Zq) is much lower.

Convergence properties with various systems of equations When mixed models are solved by iteration (GS or PCG), the convergence rate is better with fewer equations per phenotype. Experiences indicate that the convergence rate is very good with sire models, much slower with animal models, and could be especially slow with young animals not tied to phenotypes or progeny.

The regular MME equations involve an inverse of A or G, e.g.

I + A À1 α À Á â = y
An equivalent system of equation avoids creating an inverse [START_REF] Strandén | Technical note: derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit[END_REF]:

A + Iα ð Þ â = Ay
The convergence rate in such models is poor and may not occur for larger systems of equations. This is because A -1 'connects' only parents and progeny resulting in 'sparse' equations, whereas all animals in A are connected, resulting in dense equations. As G differs little from A [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Wang | Differences between genomic-based and pedigree-based relationships in a chicken population, as a function of quality control and pedigree links among individuals[END_REF], the same thinking applies to models with G and G -1 .

Computations in SNP BLUP (and BayesX)

SNP BLUP and its various forms (including BayesA, BayesB, etc.) are very popular models in genomic computations. When the number of SNP is small, the MME for SNP BLUP can be created explicitly. As the memory and compute requirements are quadratic, explicit storage of SNP BLUP equations is unfeasible with large number of SNP. For example, memory requirements for MME with 1 million SNP would be 10 12 elements also with 10 12 elements of MME created from each genotyped animal. Various options in such a case have been explored by [START_REF] Legarra | Computing strategies in genome-wide selection[END_REF] with a data set on about 2 k individuals and 20 k SNP markers. Solution by Cholesky decomposition took 2 h, and computations would increase cubically with the number of SNP, whereas memory would increase quadratically. With solution by GS and straightforward iteration on data, computations took about 4 days, with linear memory and quadratic memory costs. With a special GS implementation called GSRU, where currently computed solutions are used to adjust residuals for all individuals, computing took about 1 min and both computing and memory costs were linear. Finally, with PCG iteration on data, computing took only 20 s, with approximately linear costs for memory and computations. Although GSRU is the only realistic algorithm with millions of SNP for Bayesian models, the PCG approach can be useful for models where SNP variances are derived from SNP solutions (Zhang et al., 2010; Sun et al., 2012; Wang et al., 2012). A hybrid method would use PCG iteration Computations in genomic selection with inner (non-Monte Carlo-based) updates of weights [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF].

Variations of single-step equations

Due to the prohibitive cubic cost of inversion of G with many genotyped animals, alternative equations were proposed for ssGBLUP. Equations in [START_REF] Misztal | Computing procedures for genetic evaluation including phenotypic, full pedigree, and genomic information[END_REF] used H rather than H -1 . Convergence was achieved only with a small number of animals. Equations by [START_REF] Legarra | Computational strategies for national integration of phenotypic, genomic, and pedigree data in single-step best linear unbiased prediction[END_REF] used G and A 22 , which are easy to use, instead of G -1 and A 22 -1 . The convergence rate was slow with medium models and not achieved with large models [START_REF] Aguilar | Genetic evaluation using unsymmetric single step genomic methodology with large number of genotypes[END_REF]. Liu et al. (2014) proposed equations with SNP effects for genotyped animals. They reported lack of convergence with real data. All these models illustrate computing advantages with inverses of dispersion matrices (A -1 or G -1 ). As SNP BLUP models have good computing properties, [START_REF] Fernando | A class of Bayesian methods to combine large numbers of genotyped and non-genotyped animals for wholegenome analyses[END_REF] looked into ssGBLUP where ungenotyped animals were imputed and only SNP effects (plus imputation errors for ungenotyped animals) were estimated. Although convergence rates were good, imputations required supercomputer and imputed genotypes required large memory for big populations. Better computing times were reported using a computer with fast graphic cards [START_REF] Golden | Bolt and an alternative approach to genomic EPDs[END_REF]. [START_REF] Meuwissen | Genomic predictions based on animal models using genotype imputation on a national scale in Norwegian Red cattle[END_REF] presented a new type of H matrix based on (possibly fractional) imputation of ancestors by segregation analyses, followed by a sophisticated tailoring of G to A to compensate for errors in imputation. Although the cost of segregation analysis is high, it is not yet clear whether such an approach offers accuracy benefits when most ancestors are genotyped, and the process of tailoring is rather complex. All non-traditional ssGBLUP would require new research and investment in code for implementation of more complex models, such as multiple-trait, random regression or maternal models.

Sparse genomic relationship matrix

Computations with ssGBLUP would be easy for large number of animals if G -1 was sparse, like A -1 . Past studies suggested that G had a limited rank, as inversion is unstable with >5 to 10 k animals, and Gs as used in genetic evaluations are blended with A 22 [START_REF] Vanraden | Efficient methods to compute genomic predictions[END_REF][START_REF] Aguilar | A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF] to achieve full rank. The eigenvalue decomposition of G is:

G = UDU 0
where U is a matrix of eigenvectors and D the matrix of eigenvalues. If all eigenvalues are positive, the inverse of G is

G À1 = U 0 D À1 U
If some eigenvalues are 0, the inverse does not exist. If, due to noise, many eigenvalues are close to 0 and carry no information, they can cause instability of the inverse through large elements of D -1 . Let D t indicate a fraction of D with nonnegligible eigenvalues, and let U t be corresponding eigenvalues. Then

G À = U 0 t D À1 t U t If G -1
is to be created explicitly, the computing cost is cubic with quadratic storage. If G -is used by a PCG algorithm as

G À q = U 0 t D À1 t U t q À Á À Á
and the number of eigenvalues in D t is small, G -q can be calculated and stored at linear cost with respect to the number of genotyped animals. However, computing U and D is expensive and requires storage of G in full.

When G has a limited dimensionality, [START_REF] Misztal | Using recursion to compute the inverse of the genomic relationship matrix[END_REF] proposed to calculate G -1 using only a fraction of G and without eigenvalue decomposition. Let

G = G cc G cn G cn G nn !
where index c denotes core animals and index n non-core animals. In the algorithm for proven and young animals (APY)

G À1 % G À1 cc 0 0 0 ! + ÀG À1 cc G cn I ! M À1 ÀG nc G À1 cc I Â Ã ;
where M is a diagonal matrix with elements m i = g ii ÀG ic G À1 cc G ci for individual i in the core group. Inversion of G by APY has a cubic cost (and quadratic memory) for core individuals and a linear cost (and linear memory) for non-core individuals. Figure 1 shows an example non-zero structure of a regular inverse (after blending to make it positive-definite) and different APY inverses. The inverses are diagonal for non-core animals. As the regular G is singular for a large population, an infinite number of generalized inverses exist, including those by the APY algorithm with different sets of core animals.

Core animals can be randomly selected [START_REF] Fragomeni | Hot topic: use of genomic recursions in single-step genomic best linear unbiased predictor (BLUP) with a large number of genotypes[END_REF], with an optimal number of about 12 000 for Holsteins Misztal and Legarra [START_REF] Masuda | Implementation of genomic recursions in single-step genomic BLUP for US Holsteins with a large number of genotyped animals[END_REF] and 10 000 for Angus [START_REF] Lourenco | Genetic evaluation using single-step genomic BLUP in American Angus[END_REF]. The number of core animals is linked to effective population size (Misztal, 2016). Optimal number of core animals is equal to the number of the largest eigenvalues in D explaining 98% of the variance in G (Pocrnic et al., 2016a). With such a number, prediction of breeding values with APY was more accurate than with a regular inverse. Computation of G -1 for 570 k genotyped Holsteins by APY took only about 2 h, whereas a month computing with a large memory would be required for a regular inverse [START_REF] Masuda | Implementation of genomic recursions in single-step genomic BLUP for US Holsteins with a large number of genotyped animals[END_REF]. The APY algorithm removes computing limits from inversion of G and subsequently ssGBLUP.

Computing of A 22

ssGBLUP as in [START_REF] Aguilar | A unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score[END_REF] requires

A À1 22 . Although A -1 is sparse, A À1
22 may be relatively dense [START_REF] Faux | Inversion of a part of the numerator relationship matrix using pedigree information[END_REF], and its storage for potentially millions of genotyped animals impossible. This matrix can be calculated indirectly [START_REF] Strandén | Comparison of some equivalent equations to solve single-step GBLUP[END_REF].

A À1 22 = A 22 À A 12 À Á 0 A 11 À Á À1 A 12
where all submatrices are sparse and can be stored explicitly.

With PCG iteration, explicit A À1 22 is not required, and its product with a vector q can be calculated every round as follows [START_REF] Masuda | Implementation of genomic recursions in single-step genomic BLUP for US Holsteins with a large number of genotyped animals[END_REF]:

A À1 22 q = A 22 À A 12 À Á 0 A 11 À Á À1 A 12 h i q In particular s = A 11 À Á À1 A 12 q
is computed as a solution to:

A 11 s = A 12 q
As all matrices are sparse, the cost of above calculations is small [START_REF] Masuda | Implementation of genomic recursions in single-step genomic BLUP for US Holsteins with a large number of genotyped animals[END_REF].

GREML and APY

GREML is REML using a GRM. Similarly, REML using matrix H can be called single-step GREML (ssGREML). In the past, the success with REML in animal populations was due to sparse inverse of A -1 and sparse matrix factorization/ inversion [START_REF] Misztal | Sparse matrix inversion for restricted maximum likelihood estimation of variance components by expectationmaximization[END_REF][START_REF] Pérez-Enciso | FSPAK: an interface for public domain sparse matrix subroutines[END_REF]. If G is dense and large, [START_REF] Masuda | Technical note: acceleration of sparse operations for average-information REML analyses with supernodal methods and sparse-storage refinements[END_REF] showed that GREML with the old sparse matrix package is inefficient and unstable. A new package called YAMS recognizes dense blocks in MME, rearranges computations accordingly, allowing also for parallel computations for large dense blocks. With YAMS, GREML was much faster and more reliable.

Additional issues

This paper focused on algorithmic issues relevant to current mixed model computations in genomics in animal breeding. Many other issues exist and many will become relevant. For example, genotyping of animals with small amount of information leads to greater importance of avoidance of double counting in multi-step evaluations, or of more accurate match between genomic and additive relationships in ssGBLUP. With efforts to use sequence data, important operations are aligning short fragments, imputation, genomewide association studies, and incorporations of the detected variants in the evaluation. In dairy, a de facto global breeding scheme for Holsteins creates a need for a comprehensive genomic evaluation by Interbull that has a high accuracy, yet does not compromise privacy of the member countries. In general, while brute-force approaches may be formidable, as or more accurate yet inexpensive approaches may exist that exploit peculiarities of data including a limited effective population in farm populations. For instance, the sequence includes three G base pairs, but farm populations can be described by about 4000 to 20 000 chromosomal segments (Pocrnic et al., 2016b). This means a block size of about 150 to 750 kbases with hard to differentiate SNPs, and ability to reduce costs by working with blocks rather than individual SNP. This also suggests that increases of genomic accuracies past 4000 to 20 000 of high accuracy genotyped animals are limited (ignoring G × E and decays of predictivity over time). Similarly, if Interbull countries are unwilling to share genotypes on individual bulls or SNP effects derived from their population, they can contribute GEBV and genotypes of artificially generated animals (Maantysaari, 2015, personal communication) with their number and accuracies relevant to each country.

Conclusions

In conclusion, common operations in genetic selection are computationally feasible for an almost unlimited number of genotyped individuals using appropriate algorithms. Solving 'SNP' equations is best accomplished by iteration on data using either a special form of GS iteration or the PCG algorithm. Solving ssGBLUP is best accomplished by iteration on data using the PCG algorithm. For GREML, a sparse matrix package needs to account for dense blocks. Computations in the two latter cases greatly benefit from efficient (and sparse) inverse of the GRM.

Figure 1

 1 Figure 1 Non-zero pattern of a regular (a) and APY (b to d) inverses of the genomic relationship matrix with different choices of core animals.

https://doi.org/10.1017/S1751731116002366 Downloaded from https:/www.cambridge.org/core. INRA -Versailles-Grignon, on 01 Feb

at 10:22:04, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

https://doi.org/10.1017/S1751731116002366 Downloaded from https:/www.cambridge.org/core. INRA -Versailles-Grignon, on 01 Feb 2017 at 10:22:04, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

https://doi.org/10.1017/S1751731116002366 Downloaded from https:/www.cambridge.org/core. INRA -Versailles-Grignon, on 01 Feb 2017 at 10:22:04, subject to the Cambridge Core terms of use, available at https:/www.cambridge.org/core/terms.

Acknowledgments

Helpful suggestions from the two anonymous reviewers are gratefully acknowledged. This research was primarily supported by grants from Holstein Association USA (Brattleboro, VT), American Angus Association, Zoetis, Cobb-Vantress, Smithfield Premium Genetics, Pig Improvement Company, and the US Department of Agriculture's National Institute of Food and Agriculture (Agriculture Computations in genomic selection

and Food Research Initiative competitive grant 2015-67015-22936). The authors thank Ivan Pocrnic for help with editing.