text:
## Main

Colorectal cancer is a major health problem worldwide (Beretta et al, 2004; Bernold and Sinicrope, 2006; Monga and O’Connell, 2006). It is the second most frequent cause of cancer death in the United States and most European countries. In Japan, although its number of deaths ranks third, the disease incidence has recently been increasing (Committee of Cancer Statistics, 2005; Ministry of Health, Labor and Welfare, 2006). The prognosis of patients with colorectal cancer has improved steadily over the past several decades, owing to the development of effective therapies (Ragnhammar et al, 2001; Beretta et al, 2004; Bernold and Sinicrope, 2006; Monga and O’Connell, 2006). In practice, surgery remains the primary treatment modality for localised colorectal cancer. Adjuvant therapy is given in an attempt to eradicate micrometastases and to thereby increase the cure rate after surgical resection.

The natural history of rectal cancer differs from that of colon cancer for epidemiological, anatomical and pathophysiological reasons (Beretta et al, 2004; Bernold and Sinicrope, 2006; Monga and O’Connell, 2006). In the adjuvant therapy for stage III rectal cancer, a combined modality therapy consisting of i.v. 5-fluorouracil (5-FU) and radiotherapy is the standard care (Beretta et al, 2004; Bernold and Sinicrope, 2006; Monga and O’Connell, 2006). However, the use of chemotherapy has been based on the evidence of postoperative adjuvant chemotherapy in colon cancer, not rectal cancer (Takiuchi, 2006). In rectal cancer, there is only one adjuvant therapy (i.e., the MOF in the NSABP R-01 study), which showed improvement of overall survival (OS) and disease-free survival (DFS), compared with surgery alone as the control (Fisher et al, 1988). The optimal timing of adjuvant therapy and the combination of agents continue to evolve. As with in colon cancer (André et al, 2004; Twelves et al, 2005; Bernold and Sinicrope, 2006; Lembersky et al, 2006; Monga and O’Connell, 2006), trials of the combinations of several active drugs, such as oral fluoropyrimidines (i.e., uracil–tegafur (UFT) and capecitabine), oxaliplatin and irinotecan, are under way (Bernold and Sinicrope, 2006; Monga and O’Connell, 2006).

In Japan, in contrast to Europe and the United States above, postoperative adjuvant chemotherapy using oral fluoropyrimidines such as UFT, without radiation, has been the primary treatment, mainly because of excellent outcomes of TME with lateral pelvic lymphadenectomy (Takiuchi, 2006). The local recurrence rate in this extended TME was estimated to be lower than simple TME (Akasu and Moriya, 1997; Mori et al, 1998). Although individual RCTs comparing oral fluoropyrimidines with surgery alone showed improvement in both DFS and local recurrence in rectal cancer, but not in OS (Kodaira et al, 1998; Kato et al, 2002), meta-analyses of RCTs related to fluoropyrimidines in the 1980s showed significant benefits in both DFS and OS (Sakamoto et al, 1999; Meta-Analysis Group, 2004). In contrast, there has been no evaluation on DFS and OS in other therapies including 5-FU/leucovorin (LV) and radiation, in Japan. On the basis of the evidence, UFT is indicated as the first-line treatment of adjuvant chemotherapy after surgery for stage III rectal cancer in Japan (Japan Society for Cancer of the Colon and Rectum, 2005).

Recently, the National Surgical Adjuvant Study of Colorectal Cancer (NSAS-CC) in Japan (Akasu et al, 2006), an RCT of UFT adjuvant therapy for stage III rectal cancer, showed remarkable OS and RFS gains, compared with surgery alone. These results may be attributed to a higher dosage of UFT, rigorous surgical procedures and high quality of the clinical trial. The dosage of UFT per day of 600 mg in this study was greater than the 300–400 mg used in previous UFT studies (Kodaira et al, 1998; Sakamoto et al, 1999; Kato et al, 2002; Meta-Analysis Group, 2004). In this study, TME with lateral pelvic lymphadenectomy was chosen as the control treatment.

Under these circumstances of treatment evolution in rectal cancer, patients and physicians have to consider risks and benefits in choosing the best treatment between the different options now available. On the other hand, health-care payers increasingly require evidence of economics in addition to the clinical value, under the severe pressure to health-care expenditure. However, there is far less economic evaluation specifically of adjuvant therapy for rectal cancer than for colon cancer (van den Hout et al, 2002). Only one study, related to cost-effectiveness of adjuvant therapy for rectal cancer, has been reported (Norum et al, 1997). Our objective was to confirm the cost-effectiveness of adjuvant UFT therapy in Japan, extrapolating these results to a time horizon sufficiently distant to capture all costs and outcomes of relevance. This study would provide basic information on the cost-effectiveness of adjuvant therapy for rectal cancer not only in Japan, but also in other countries indirectly.

## Materials and methods

### Analytical overview

Economic analysis was conducted retrospectively based on the Japanese NSAS-CC (Akasu et al, 2006), a multicentre RCT. Patients with completely resected stage III rectal cancer, who underwent standardised TME with selective lateral pelvic lymphadenectomy, were randomly assigned to either oral UFT (400 mg m−2 per day) for 1 year (n=139) or surgery alone (n=135).

As a type of economic analysis (Drummond et al, 2005), a cost-effective analysis was performed. Incremental costs and effectiveness of adjuvant UFT therapy compared with surgery alone were evaluated. According to the effectiveness measure used (i.e., life-years gained and quality-adjusted life-years (QALYs) gained), incremental cost-effectiveness ratios (ICERs) were calculated based on cost per life-year gained and cost per QALY gained, respectively. In a narrow meaning, the former is termed cost-effectiveness analysis, whereas the latter is termed cost–utility analysis.

The payer of National Health Insurance in Japan was adopted as a perspective of economic analysis (Drummond et al, 2005). Therefore, for costs, direct medical-care costs (e.g., costs of tests, drugs, health-care personnel and so on) were examined, whereas indirect costs (e.g., time costs or production loss among patients and their families) were not considered. As a time horizon for evaluation, three levels of time periods (i.e., observational period (5.6 years), 10-year follow-up and over lifetime) were considered. As the base case analysis, observational period was used, as this period incorporated few assumptions for the evaluation.

### Effectiveness

The results of the NSAS-CC were used as evidence of effectiveness in the economic analysis. The clinical results have been presented in detail elsewhere (Akasu et al, 2006). As is shown in Table 1, between the UFT therapy group and the surgery alone group, no statistical differences were observed in age, sex, tumour location or pathological tumour stage. The incidence of adverse events more than grade 3 in the UFT therapy group was higher than that in the surgery alone group. The OS and RFS rates in the UFT therapy group were higher than those in the surgery alone group (Akasu et al, 2006).

Using patients’ data, OS and RFS were estimated by the Kaplan–Meier method, up to 5.6 years from randomisation. Beyond the observation period of 5.6 years, OS was simulated using the Boag model (Boag, 1949) combined with the competing risk model (Gross and Clark, 1975; Maetani et al, 2004) (Figure 1). Commonly used methods for extrapolating the survival data beyond the observation include Weibull, Gompertz, exponential, log-normal and generalised gamma distributions (Lee and Go, 1997). However, there is no explicit standard for deciding the optimal method. In this analysis, the Boag model combined with the competing risk model was used. This model consisted of two components: the disease-specific survival curve and the disease-independent survival curve. Although the former was simulated by a log-logistic model, the latter was simulated by the survival curve of the general population matched for age and sex of the subjects. Instead of the log-normal model adopted in the original Boag model, the log-logistic model was selected in this analysis, according to Akaike's Information Criteria (Akaike, 1974).

The mean number of life-years and relapse-free life-years for patients in each group was estimated as the area under RFS and OS curves (Willan and Briggs, 2006). In addition, QALYs were estimated from OS and RFS by weighting each survival by utility value for each possible health state. Utility values for the health states were derived from the published study by Ness et al (1999). This study was identified by a systematic literature search using MEDLINE and EMBASE and critically appraised. The median utility (25th, 75th percentile) for stage III after surgery and metastasis was 0.75 (0.55, 0.85) and 0.20 (0.00, 0.40), respectively. These utility values were assigned to the health conditions with and without relapse in this analysis. Although the frequency and grade of adverse events due to chemotherapy was relatively low, utility reduction associated with these events was adjusted by the method adopted by Aballéa et al (2007).

### Cost

Costs incurred for resources used during trial and subsequent follow-up were estimated from trial data and their extrapolation. Resource utilisation during trial and follow-up was derived from individual patient history data. As observations on many patients are censored in a clinical trial, subsequent costs are unknown. To correct for censoring, the cost history method proposed by Lin et al (Lin et al, 1997; Willan and Briggs, 2006) was applied. Costs were estimated from the National Health Insurance perspective using the National Health Insurance reimbursement list and drug price in 2005 (Institute of Social Insurance, 2005; Jiho, 2005). Unit costs for chemotherapy (UFT), consultation and diagnostic tests were low, whereas those of imaging tests for follow-up were relatively high. Admission fees were considerably lower than those in the United States.

### Sensitivity analysis

The results of a probabilistic sensitivity analysis are shown in Figure 2. Figure 2A shows ICER (cost per QALY gained) scatter plots. More than 95% of the points resided in the southeast quadrant (i.e., more effective and less costly). The cost-effectiveness acceptability cure based on 5000 samples of cost-effectiveness ratio is presented in Figure 2B. Even if additional QALY was valued as 0, the likelihood of UFT therapy being cost-effective was 98%. The NMB curve is shown in Figure 2C. The NMB and its confidence interval (CI) curves did not cross the horizontal axis. This indicates that UFT therapy was beneficial, even if a decision maker was not willing to pay anything for the additional QALY.

A number of qualitative sensitivity analyses are shown in Tables 2B and 4. As to time horizon (Table 2B), from 5.6-year observation to over lifetime, cost-effectiveness ratios were all negative, indicating more benefits and less costs (i.e., dominance).

The two-way sensitivity analysis of discount rate for both costs and effect did not show any change in the dominance of UFT therapy. Incremental cost-effectiveness ratio was lowest (−$5031 per QALY) at a discount rate of 5% for both costs and effectiveness and highest (−$4933 per QALY) without discounting. Incremental cost-effectiveness ratio increased with increase of discount rate of cost, whereas ICER decreased with increase of discount rate of effect.

The results of one-way sensitivity analyses are shown in Table 4. Variations in recurrence rate, utility value, QALYs and the acquisition cost of UFT did not influence the dominance of UFT therapy. On the other hand, with variations of total cost, ICERs varied from −$9695 to$334 per QALY gained. This upper limit of ICER was positive, but at a very low level.

### Budget impact

The annual incidence rate of rectal cancer was 3.23 × 104 and the proportion of the stage III patients was 32.5%. Then, the annual number (95% CI) of the resected stage III rectal cancer was estimated to be 1.05 × 104 (0.91–1.20). The cost reduction per patient was $2603. Therefore, a budget impact converting from surgery alone to adjuvant UFT therapy would reduce current medical expenditure by$27.3 million (95% CI, $23.7 million,$31.2 million) in the first 5.6 years after its adoption.

## Discussion

From the perspective of the National Health Insurance in Japan, this cost-effectiveness analysis of UFT adjuvant therapy for stage III rectal cancer would save health-care costs and improve health outcomes, compared with surgery alone (Table 2). Uracil–tegafur therapy has proved dominant (less costly and more effective). The cost–utility ratio of UFT can be ranked near the top of league table of cost–utility in oncology (Earle et al, 2000). The cost savings on the National Health Insurance budget during the observational period would be approximately $27.3 million. There has been little evidence on economic evaluation of adjuvant therapy for rectal cancer (van den Hout et al, 2002). Only one cost-effective analysis of adjuvant therapy for Dukes’ B and C colorectal cancer (Norum et al, 1997), of which about two-third consisted of colon cancer, has shown that cost per QALY gained of 5-FU/levamisole was between £4800 and £16,800, compared with surgery alone. This result is highly cost-effective compared with the recent threshold (Jonsson, 2004) for cost-effectiveness. However, cost-effectiveness ratio in our study is much better than that, even though these results are not directly comparable with each other. In the area of rectal cancer, economic evaluations on the standard adjuvant therapy (i.v. 5-FU/radiation) comparing control (i.e., surgery alone), as well as head-to-head evaluation comparing newly emerging with the standard therapy, are urgently needed for health-care decision making. Our study provides basic information in comparing cost-effectiveness of adjuvant therapies for rectal cancer. To estimate stochastic uncertainty of ICER due to sampling bias in this study, probabilistic sensitivity analyses (Glick et al, 2001; Briggs, 2004; Drummond et al, 2005) were performed (Table 2, Figure 2). Cost-effectiveness scatter plots showed that CIs of ICER located in the southeast quadrant (more effective and less costly) on the cost-effectiveness plane. Cost-effectiveness acceptability and NMB curves give more information than simple intervals or plots of ICERs mentioned above. Even if a decision maker was unwilling to invest anything at the maximum to achieve additional QALY, the likelihood of UFT therapy being acceptable as cost-effective was 98% (Figure 2B). At the same condition, the NMB curve showed that UFT therapy was beneficial ($2450 per person) (Figure 2C). These results show that the dominant cost-effectiveness of UFT adjuvant therapy is robust. This dominance was shown more sensitive to costs than effectiveness by the scatter plot of ICERs (Figure 2A) and broad one-way sensitivity analysis on both QALY gained and total costs (Table 4).

The time horizon is an important issue to sufficiently capture relevant costs and health outcomes of UFT adjuvant therapy. The observation period of the NSAS-CC, 5.6 years, was limited. Although most costs were incurred mainly in the observational period, life-years gained would continue after it. In this study, a simulation model was used to extrapolate its results. There is a variety of ways for simulation (Lee and Go, 1997), but no uniform methodology is available. In the analysis, we used the Boag model combined with the competing risk model (Maetani et al, 2004). In a sensitivity analysis, ICER of the observational period (the base case analysis) was compared with that in the extrapolated periods (i.e., 10-year follow-up and over lifetime). These ICERs were dominant (more effective and less costly) in spite of relative differences in their values.

The key drivers of the dominant cost-effectiveness results of UFT are mainly the savings achieved by reduction of costs related to recurrence and death, which offset the acquisition cost of UFT. After a 4-year follow-up, any recurrence has not been observed in this study. In one-way sensitivity analysis (Table 4), varying recurrence rates between 95% CIs did not have any substantial impact on dominant cost-effectiveness. The other cost driver was the acquisition costs of UFT. Variation in the price of UFT and its standard regimen would be unlikely. Varying acquisition costs of UFT did not have any impact of dominance (Table 4).

Cost-effectiveness analysis using QALYs offers the opportunity to consider both quantity and quality of survival. However, no substantial difference in ICERs was observed between cost per LY gained and QALY gained (Table 2). In addition, the sensitivity analysis on range of utility values for recurrence and non-recurrence revealed no major change in dominant cost-effectiveness (Table 4). In a sizable fraction of cost-effectiveness analysis, utility weighting was indicated not to substantially alter the estimated cost-effectiveness of an intervention (Chapman et al, 2004). It is thus suggested that sensitivity analyses using ad hoc adjustment or weight from the literature may be sufficient. Our results support this conclusion.

The impact of discounting for the time value of money on the results was examined extensively by two-way sensitivity analyses. Although ICERs were more sensitive to cost discounting than effectiveness discounting, there was no substantial change in dominant cost-effectiveness. The main reason is likely to be that major costs were incurred during the early phase of follow-up and improved survival was realised simultaneously.

There are several limitations in the analysis that should be commented on, and the results should be treated with caution. First, the analysis was based on a small RCT with relatively short-term follow-up (Akasu et al, 2006). Therefore, extensive sensitivity analyses were performed to examine this uncertainty. However, large-scale RCTs are crucial to resolving this issue. Second, the UFT adjuvant therapy was the only chemotherapy that showed improvement in OS and DFS compared with surgery alone in Japan. As there is no RCT comparing the standard adjuvant therapy (e.g., 5-FU/radiation) with surgery alone or UFT, it is impossible to directly or indirectly compare these therapies. In future, head-to-head evaluations comparing new emerging therapies (e.g., capecitabine and oxaliplatin) will need to be carried out. Third, the perspective of this analysis is that of a payer for health care, rather than a society. From a societal perspective, the range of costs is broader and includes time costs and travel costs associated with treatment and loss of production due to earlier death. As UFT adjuvant therapy increased OS and decreased recurrence, their cost reduction in both indirect costs and direct costs will offset the costs corresponding to a treatment period.

The issue of generalisability of this study to other countries should be carefully examined, as the UFT adjuvant therapy is the standard in Japan. However, as mentioned in the introduction, in contrast to colon cancer, there is no firm evidence for effectiveness of the western standard therapy (e.g., FU/LV) in rectal cancer. Moreover, given the high treatment costs, substantial toxicity and relatively limited efficacy of the fast-changing chemo- and immunotherapeutic combinations for colorectal cancer, examination of cost-effectiveness studies should be conducted on a routine basis along with determination of clinical benefits (Jansman et al, 2007). Uracil–tegafur has been approved and utilised in 31 countries from Europe and Canada to Asia, excluding the United States. Recently, several studies of combination therapies of UFT/LV in metastatic colorectal cancer have shown efficacy (Bennouna et al, 2006; Bajetta et al, 2007). Therefore, the results of this study indicating effectiveness and cost-effectiveness of UFT adjuvant therapy in rectal cancer may be useful as a reference case for the direct or indirect future examination internationally.


text:
On Linear Convergence of Weighted Kernel Herding

# On Linear Convergence of Weighted Kernel Herding

Rajiv Khanna Department of Statistics, University of California at Berkeley, Email: rajivak@berkeley.edu      Michael W. Mahoney ICSI and Department of Statistics, University of California at Berkeley, Email: mmahoney@stat.berkeley.edu
###### Abstract

We provide a novel convergence analysis of two popular sampling algorithms, Weighted Kernel Herding and Sequential Bayesian Quadrature, that are used to approximate the expectation of a function under a distribution. Existing theoretical analysis was insufficient to explain the empirical successes of these algorithms. We improve upon existing convergence rates to show that, under mild assumptions, these algorithms converge linearly. To this end, we also suggest a simplifying assumption that is true for most cases in finite dimensions, and that acts as a sufficient condition for linear convergence to hold in the much harder case of infinite dimensions. When this condition is not satisfied, we provide a weaker convergence guarantee. Our analysis also yields a new distributed algorithm for large-scale computation that we prove converges linearly under the same assumptions. Finally, we provide an empirical evaluation to test the proposed algorithm for a real world application.

## 1 Introduction

Estimating expectations of functions is a common problem that is fundamental to many applications in machine learning, e.g., computation of marginals, prediction after marginalization of latent variables, calculating the risk, bounding the generalization error, estimating sufficient statistics, etc. In all these cases, the goal is to approximate tractably the integral

 Ep[f(x)]=∫Xf(x)p(x)dx, (1)

for a given function , with value oracle access, and density , with sampling oracle access. In most cases of interest, the integral is not computable in closed form. A common solution is to sample from the domain of , and then to approximate the full integration using the samples. The simplest version of this is to use vanilla Monte Carlo (MC) Integration—sample uniformly at random from the domain of , and then output the empirical average as the estimate of the integral. This simple algorithm converges at a prohibitively slow rate of . More generally, one can apply Markov Chain Monte Carlo (MCMC) algorithms for (1), and these algorithms come with many benefits, but they also converge quite slowly.

To alleviate the slower convergence of MC-based Integration, non-uniformly at random sampling algorithms such has Kernel Herding (KH) have been proposed. The initial herding algorithm [18, 19, 20] was proposed to learn Markov Random Fields (MRFs), and it was applicable to discrete finite-dimensional spaces. This was extended to continuous spaces and infinite dimensions using the kernel trick by [3], who also provided a convergence rate of , which is faster than the convergence rate of MC Integration. The speed up can be attributed to a smarter selection of sampling points (rather than uniformly at random) that ensures that the selected samples are somewhat negatively correlated. In addition, [3] also provided a moment matching interpretation of the algorithm. Essentially, solving (1) using the weighted average of a few samples () can be interpreted as minimizing the discrepancy in a Reproducing Kernel Hilbert Space (RKHS). The discrepancy is also widely known to be the Maximum Mean Discrepancy (MMD) metric [3, 8].

In an alternative Bayesian approach, [15, 6] assume a Gaussian Process prior on and use a herding algorithm to sample so as to minimize the corresponding moment matching discrepancy arising out of the posterior of , after observing the points already sampled. This perspective was shown to be equivalent to minimizing the herding discrepancy, but with an additional step of also minimizing the weights attached to the sampled points, instead of using uniform weights of  [8]. Empirically, it was observed that this new algorithm—called Sequential Bayesian Quadrature (SBQ)—converges faster than naive herding because of the additional weight optimization step. To the best of our knowledge, the analysis for faster convergence for SBQ was still largely an open problem before our present work.

A partial explanation for faster convergence in Weighted Kernel Herding (WKH) was provided by [2]. While SBQ chooses the next sample point and the weights to minimize the said discrepancy, the herding algorithm only solves a linearized approximation of the same discrepancy. Once, the sample point is selected, the weights themselves can be optimized. We name this algorithm (solving linearized discrepancy + optimizing weights) as WKH. For the case when the weights are constrained to sum up to [2] noted that WKH is equivalent to the classic Frank-Wolfe [5] algorithm on the marginal polytope. Exploiting this connection, they were able to use convergence results from Frank-Wolfe analysis to analyze constrained WKH. Specifically, if the optimum lies in the relative interior of the marginal polytope, so that it is distance away from the boundary, the convergence rate is . For finite dimensional kernels, they show ; while for infinite dimensional kernels, . Even for finite dimensions, can be arbitrarily close to , rendering the exponential convergence meaningless. It was also pointed out by [2] that the existing theory does not fully justify the much better empirical performance of WKH. To the best of our knowledge, this question has remain unresolved before our present work.

In our present work, we provide an analysis for exponential convergence of unconstrained WKH. As noted by [8], the weights in WKH do not need to sum to in many applications, and so the correspondence to Frank-Wolfe and its limitations do not hold for us. Furthermore, our analysis also extends to non-convex support of and infinite dimensional kernels, provided a simplifying structural assumption is satisfied.

Central to our analysis is what we call a realizability assumption. We assume that the mean in (1) can be exactly reproduced by a linear combination of samples , from the domain. We discuss the assumption and its impact in Section 3, and we emphasize that it is not restrictive. In fact, for many nice cases, , but more generally, realizability guarantees convergence. The comparative view of the algorithms and their respective rates are presented in Table 1.

A side effect of our analysis is that we are able to scale up the herding algorithm to multiple machines, cutting down on run time and memory requirements. As noted by [8], the iteration of KH and WKH require an search over the domain to sample the next point, while SBQ takes . Thus, for larger , these algorithms can get slower. Our proposed solution is to split the domain and hence the search to multiple machines, and to run a local WKH/SBQ, thus speeding up the search for the next sample point. The local solutions are then collated to output the final set of samples. Interestingly, we show that this distributed algorithm also exhibits linear convergence.

Our contributions in this work are: (a) Theoretical—providing the fastest known convergence rates for two popular sampling algorithms; (b) Algorithmic—providing a new distributed herding algorithm that we also show converges linearly; and (c) Empirical—since there is ample empirical evidence supporting the good performance of the KH/SBQ algorithms in earlier works [8, 3, 2], we focus on and demonstrate the empirical performance of the distributed algorithm in a real world application. In more detail, our main contributions are:

• We analyze and prove linear convergence of two widely used algorithms for approximating expectations, under standard assumptions for well-known setups. Our analysis also provides theoretical justifications for empirical observations already made in the literature.

• We propose the novel niceness assumption of realizability in the optimization problem we study. The assumption trivially holds for many known cases, but even for several harder cases, it acts as a sufficient condition to guarantee linear convergence.

• We consider pathological infinite dimensional cases, and we provide a weaker convergence guarantee for cases that do not follow realizability.

• We propose a novel distributed algorithm for approximating expectations for large scale computations, and we study its convergence properties.

• Finally, we present empirical studies to validate the newly proposed algorithm.

#### Other related works.

The Frank-Wolfe algorithm [10, 11] has other variants than the one considered by [2] to relate to KH. These variants enjoy faster convergence at the cost of additional memory requirement. Specifically, instead of just selecting sample points, one can think of removing bad points from the set already selected. This variant of Frank-Wolfe, known as FW with away steps, is one of the more commonly used ones in practice, because it is known to converge faster. If the weights are not restricted to lie on the simplex, the analogy to matching pursuit algorithms is obvious. We refer to [14] for corresponding convergence rates. However, the linear rates for these algorithms require bounding certain geometric properties of the constraint set, and this may not be easy to do for RKHS-based applications that usually employ KH and/or SBQ. With the goal of interpreting blackbox models, [16] recently exploited connections with submodular optimization to provide the weaker type convergence we discuss (in Section 3.2) in the form of an approximation guarantee for SBQ for discrete . Our result is more general, since we also address the case of discrete (in Section 3.2). We do not make use of submodular optimization results, but instead we directly employ simple linear algebra after the key insight of realizability. A similar proof technique was also used by [17] for proving approximation guarantees of low rank optimization. The proof idea for the distributed algorithm was inspired from tracking the optimum set, which is a common theme in analysis of distributed algorithms (see, e.g., [1] and references therein).

#### Outline.

The rest of the paper is as follows. We present relevant background on the algorithms in Section 2. The main theoretical results are presented in Section 3, albeit most of the proofs are relegated to the appendix; and the algorithm for distributed kernel herding is presented in Section 4. Finally, we provide empirical results in Section 5; and we provide a brief conclusion in Section 6.

## 2 Background

In this section, we discuss some relevant background for the algorithms and methods at hand. We begin by establishing some notation. We represent vectors as small letter bolds, e.g., . Matrices are represented by capital bolds, e.g.. . Matrix transposes are represented by superscript . We use to represent probability densities over random variables, which may be scalar, vector, or matrix valued, which shall be clear from context. Sets are represented by sans serif fonts, e.g., ; and the complement of a set is . A dot product in an RKHS with kernel is represented as , and the corresponding norm is . The dual norm is written as . We denote as .

#### Weighted Kernel Herding.

Our goal is to evaluate the expectation of a function under some distribution . In general, such evaluations may not be tractable if we do not have access to the parametric form of the said density and/or function. In this work, we only assume value oracle access to the function and a sampling oracle access to the distribution. Thus, we intend to approximate the expectation of a function by a weighted sum of a few evaluations of the said function. Say a function is defined on a measurable space . Consider the integral:

 Ep[f(x)]=∫Xf(x)p(x)dx≈n∑i=1wif(xi), (2)

where are the weights associated with function evaluations at . The algorithms used to solve (2) have a common underlying theme of two alternating steps: (1) Generate the next sample point ; and (2) Evaluate the weights for all the samples chosen so far. For example, using and sampling uniformly at random over recovers the standard MC integration. Other methods include Kernel Herding (KH [3]) and quasi-Monte carlo [4], both of which use but use specific schemes to draw . Sequential Bayesian Quadrature (SBQ [15]) goes one step further and employs a non-uniform weight vector as well. For convenience, define the set . Say .

We present a brief overview of these algorithms next, and point the reader to [8, 3] for further details. As already noted, solving (2) can be equivalently considered to be minimizing a discrepancy in an RKHS. Kernel Herding (KH) chooses the next sample by solving the linearized form of the objective discrepancy, and it sets all the weights to be . Say is the dot product in the RKHS defined by the kernel function . If is the current estimate of the integral (2), then the linear objective to solve for the next sample point is given by:

 xKHn+1←argminx∈X⟨Z(Sn),x⟩k=argminx∈Xn∑i=1wik(x,xi). (3)

The SBQ algorithm is more complicated. It assumes a functional Gaussian Process (GP) prior on with the kernel function . This implies that the approximate evaluation in (2) is a random variable. The sample is then chosen as the one that minimizes the posterior variance while the corresponding weights are calculated from the posterior mean of this variable. The details follow.

Say we have already chosen points: . From standard properties of Gaussian processes, the posterior of , given the evaluations , has the mean function:

 ^f(x)=k⊤K−1f,

where is the vector of function evaluations , is the vector of kernel evaluations , and is the kernel matrix with . The quadrature estimate provides not only the mean, but the full distribution as its posterior. The posterior variance can be written as:

 cov(x,y)=k(x,y)−k(x,X)K−1k(X,y),

where is the matrix formed by stacking , and the kernel function notation is overloaded so that represents the column vector obtained by stacking . The posterior over the function also yields a posterior over the expectation over defined in (2). Then, it is straightforward to see , where . Note that the weights in (2) can be written as . We can write the variance of as:

 var(Z(Sn))=∬k(x,y)p(x)p(y)dxdy−z⊤K−1z. (4)

The algorithm Sequential Bayesian Quadrature (SBQ) samples for the points in a greedy fashion with the goal of minimizing the posterior variance of the computed approximate integral:

 xSBQn+1←argminx∈Xvar(Z(Sn∪{x})). (5)

It has been shown that SBQ also can be seen to be choosing samples so as to minimize the same discrepancy as the one by KH. Instead of solving a linearized objective as done by KH, it chooses the next sample so as to directly minimize the discrepancy i.e. .

In this work, we first analyze Weighted Kernel Herding (WKH) (Algorithm 1) which combines the best of the two algorithms of SBQ and KH by using , with updating weights using SBQ’s . The linearized objective was also considered by [2], with the additional constraint that the weights are positive and sum to . They use the connections to the classic Frank-Wolfe algorithm [5] to establish convergence rates under this additional constraint. However, their linear convergence rate is dependent on how far inwards does the optimum lie in the marginal polytope, and it devolves to the sublinear rate when the optimum can be arbitrarily close to the boundary of the marginal polytope. As noted by [8], the additional condition on the weights is not required for many practical applications. We must note, however, that in some applications it is imperative that the iterates lie within the marginal polytope [13, 12].

## 3 Convergence Results

In this section, we establish linear convergence of WKH. The starting point of our analysis is the following re-interpretation of the posterior variance minimization (4) as a variational optimization of Mean Maximum Discrepancy (MMD) [8]. We can re-write (4) as a function of a set of chosen samples as:

 \vspace−4mmg(S):=∥μp−∑i∈Swiϕ(xi)∥2k, (6)

where is the respective feature mapping, i.e., , and where is the mean embedding in the kernel space.

Before presenting our result, we delineate the assumptions we make on the cost function (6).

#### Assumption 1 (Realizability):

We assume a set of atoms for , such that .

Assumption 1 posits there exists a set of samples in the domain whose weighted average exactly evaluates to the expectation (2) of under . We emphasize that this assumption is not restrictive. Consider the case for finite dimensional , convex domain of continuous density , and continuous . Then, the expectation will lie within . Hence, just from continuity properties, the average (2) must be obtained somewhere on the domain for some , and so . For non-convex , the expectation might lie outside , but a linear combination of small number of points in may still reproduce this value. For discrete , and/or for infinite dimensional , if Assumption 1 is satisfied, we prove that linear convergence holds. To the best of our knowledge, such a sufficient condition for these harder cases has not been established before, and this is one of our contributions of this work.

We assume that the eigenvalues of the kernel matrix formed by applying the kernel function on are lower bounded by and upped bounded by .

Assumption 2 is a standard assumption of restricted strong convexity (RSC) and smoothness (RSS) over the cost function (6).

#### Assumption 3 (Standardization):

We assume that the feature mapping is standardized, i.e., , as assumed in previous works [8] for ease of exposition.

###### Theorem 1.

Under Assumptions 1 through 3, if is the sequence of iterates produced by Algorithm 1, the function converges as .

Proof Outline. The central idea of the proof is to track and bound the selection of a sample at each iteration, compared to the ideal selection that could have provided optimum solution. For this purpose, the properties of the selection subproblem (3) and the assumptions are used. The detailed proof is presented in the appendix.

### 3.1 Relationship with SBQ

The SBQ Algorithm [6, 8] is similar to the WKH algorithm, except that it solves a slightly more involved optimization subproblem to select the next sample point. The algorithm pseudocode is exactly the one presented in Algorithm 1, with step 4 replaced with (5) instead of (3). Thus, WKH solves a linear program every iteration, while SBQ solves a kernel minimization problem. In this and other theoretical works, it is assumed that the oracle that solves the optimization subproblem is given and the convergence rates generally consider the number of calls made to this oracle for desired accuracy of the overall solution. In application, however, the subproblems can be vastly impactful on convergence, run times, and memory management. It has been observed that SBQ converges faster in practice, which is intuitive since it performs more work per iteration. Furthermore, for distributed algorithm variants like the one presented in Section 4, the cost per iteration may not be a bottleneck due to availability of multiple machines. On the other hand, in addition to higher cost per iteration, the application of SBQ might not be always feasible. We present one such example in Section 5, where SBQ requires memorization of the kernel matrix for a reasonable run time.

While the two algorithms have their specific advantages, the decrease in the cost function (6) per iteration of SBQ is more compared to its decrease per iteration of WKH. This relationship, combined with Theorem 1 establishes linear convergence of SBQ, which to the best of our knowledge, was still an open question (see Table 1 in [8] which lists the convergence rate of SBQ as unknown).

###### Corollary 2.

Under Assumptions 1 through 3, if is the sequence of iterates produced by SBQ, the function converges as .

###### Proof.

The proof follows from the fact that SBQ converges at least as fast as WKH, since SBQ makes faster progress per iteration. To see this, consider a single iteration of each from a given set of iterates . Let and be the corresponding sample selected through an iteration of SBQ and WKH respectively. Then, from (5),

### 3.2 Breaking realizability

While realizability trivially holds with for many nice cases, it is important to discuss cases when realizability breaks. We first consider the case of discrete . Say is the size of the support of discrete , and is the dimensionality of the space. If , from the fact that each iteration only adds a sample that is in the orthogonal complement of previously added samples, . In this case, the same rate will still hold. On the other hand, it is possible that the mean is realizable for no less than samples if . In this case, the linear rate is no longer valid. One can call upon the Frank-Wolfe connection [2], since as long as , we know that the a linear rate holds (because the optimum lies in the relative interior away from the boundary, since the effective rate is minimum of the two). However, as discussed before, the optimum may lie arbitrarily close to the boundary. In such cases, the following statement still holds nonetheless. Say , then for iterations of Algorithms 1 or 2 (WKH or SBQ), to select the will have .

Things are trickier for infinite dimensional kernels. It is possible that in infinite dimensions, is not realizable for any . We characterize this as a pathological case, since it implies that is not in the span of any finite number of atoms . From Theorem 1, one can, however, provide the following weaker guarantee. For any projection of onto a finite number of atoms , we need to run atmost iterations of Algorithms 1 or 2 with WKH or SBQ, in order to get close to the objective value of the best such -sized projection. The proof of this statement follows from replacing with the projection in the analysis of Theorem 1.

## 4 Distributed Kernel Herding

In many cases, the search over the domain can be a bottleneck causing the herding algorithm to be too slow to be useful in practice. In this section, we develop a new herding algorithm that can be distributed over multiple machines and run in a streaming map-reduce fashion. We also provide convergence analysis for the same, using techniques developed in Section 3.

The algorithm proceeds as follows. The domain is split onto machines uniformly at random. Each of the machines has access to only , such that and for . Each machine runs its own herding algorithm (Algorithm 1) independently, by specializing (3) to do a restricted search over , instead of . Finally, the iterates being generated by each machine are sent to a central server machine which collates the samples from different machines by running another copy of the same algorithm, with replaced by the discrete set of samples sent to the central server by the other machines. Finally, the best solution out of all the solutions is returned. The pseudo-code is illustrated in Algorithm 2.We now provide the convergence guarantees for this algorithm.

###### Theorem 3.

Under Assumption 1, with , and Assumptions 2 and 3, if is the sequence of iterates produced by Algorithm 2, the function converges as .

Proof Idea: Note that the final set of filtered iterates outputted are the best out of possible sequences. The proof tracks the possibilities for when is split. The goal is to then show that under all possible scenarios, at least one of the sequences converges linearly. The convergence of individual sequences is based on the proof techniques used in proof of Theorem 1.

We remark that, for the more general case of , our proof technique does not give a valuable convergence rate. In particular, our analysis yields , which is not very meaningful. This could be an artifact of our proof technique, and improving it is an interesting open question for future research. Also, while we presented the analysis for Distributed WKH, by Corollary 2, the same rate also holds for Distributed SBQ.

## 5 Experiments

We refer the readers to earlier works for empirical evidence of the good performance of the Kernel Herding and SBQ algorithms [2, 3, 20, 18, 19, 8, 6]. In this section, we focus on studying these algorithms and the distributed variant on a real world problem of data summarization using coresets on three different datasets of varying sizes.

The task of data summarization is as follows. Our goal is to select a few data samples that represent the data distribution sufficiently well, so that a model built on the selected subsample of the training data does not degrade too much in performance on the unseen test data. More specifically, we are interested in approximating the test distribution (i.e., discrete ) using a few samples from the training set. Hence, algorithms such as SBQ and WKH are applicable, provided we have a reasonable kernel function. We use Fisher kernels [9], since they were recently used to show strong performance for this task [16]. For completeness, we provide a brief overview of the Fisher kernel in the appendix.

Another method that also aims to do training data summarization is that of coreset selection [7], albeit with a different goal of reducing the training data size for optimization speedup while still maintaining guaranteed approximation to the training likelihood. Since the goal itself is optimization speedup, coreset selection algorithms typically employ fast methods, while still trying to capture the data distribution by proxy of the training likelihood. Moreover, the coreset selection algorithm is usually closely tied with the respective model, as opposed to being a model-agnostic.

We employ different variants of WKH/SBQ to the problem of training data summarization under logistic regression, as considered by [7] using coreset construction. We experiment using three datasets ChemReact, CovType and WebSpam. ChemReact consists of chemicals each of feature size . Out of these, are test data points. The prediction variable is and signifies if a chemical is reactive. CovType has examples each of feature size . Out of these, are test points. The task is to predict whether a type of tree is present in each location or not. WebSpam has 350,000 webpages each having features. Out of these, are test data points. The task here is to predict whether a webpage is spam or not. We refer to [7] for source of the datasets.

In each of the datasets, we further randomly split the training data into validation and training. We train the logistic regression model on the new training data, and we use the validation set as a proxy to the unseen test set distribution. We build the kernel matrix and the affinity vector , and we run different variants of sampling algorithms to choose samples from the training set to approximate the discrete validation set distribution in the Fisher kernel space. Once the training set samples are extracted, we rebuild the logistic regression model only on the selected samples, and we report negative test likelihood on unseen test data to show how well has the respective algorithm built a model specific dataset summary.

The algorithms we run are WKH, SBQ, -SBQ and -WKH, where represents the number of machines used to select samples for different values of . The smaller ChemReact data fitted on a single machine, so we run WKH and SBQ on single machine. To present the tradeoff, we also run -WKH and -SBQ. These are about times faster than their single machine counterparts, but they degrade in predictive performance. For the baselines, we use the coreset selection algorithm and random data selection as implemented by [7]. The results are presented in Figure 1. We note that our algorithm yields a significantly better predictive performance compared to random subsets and coresets [7] with the same size of the training subset across different subset sizes. We note that generally SBQ has better performance numbers than WKH for same across different values of . Note that WebSpam and CovType were too big to run on a single machine, and they are thus perfect examples to illustrate the impact and usefulness of the distributed algorithm. The m-WKH-1 experiment presented in Figure 1 is obtained by averaging the output of a single split of the dataset run on one of the machines (as was done by [16] to scale up sine they do not have the distributed algorithm to work with). All the experiments were run on 12-core 16Gb RAM machines.

An interesting practical consideration between SBQ and WKH is that the inner subproblem (5) is harder than (3), not just in time complexity but also in its memory requirement. This is because solving  (5) will make repeated calls to calculate dot products across iterations, and this can rack up the run time quickly if the kernel matrix is pre-computed and stored in memory. On the other hand, (3) does not have this shortcoming, and the kernel dot product can be calculated as required without significantly impacting the runtimes.

## 6 Conclusion

We have presented novel analysis for two commonly used algorithms as well as for a new distributed algorithm for estimating means. Our results bridge the gap between theory and empirical evidence by presenting the fastest known convergence rates for these algorithms. Our realizability assumption is the key insight that allows us to improve upon previous results, and at the same time we provide sufficient conditions for convergence in the pathological infinite dimensional case. There are quite a few interesting future research directions. An important one is studying convergence for the distributed algorithm for the case with realizability . Can the realizability help in proposing other faster variants? Alternatively, what are the best possible information theoretic rates for these problems? We hope our work inspires future works in these and other related directions.

## References

• [2] Francis Bach, Simon Lacoste-Julien, and Guillaume Obozinski. On the equivalence between herding and conditional gradient algorithms. In Proceedings of the 29th International Coference on International Conference on Machine Learning, ICML’12, pages 1355–1362, 2012.
• [3] Yutian Chen, Max Welling, and Alexander J. Smola. Super-samples from kernel herding. In UAI, 2010.
• [4] Josef Dick and Friedrich Pillichshammer. Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration. Cambridge University Press, New York, NY, USA, 2010.
• [5] Marguerite Frank and Philip Wolfe. An algorithm for quadratic programming. Naval research logistics quarterly, 1956.
• [6] Zoubin Ghahramani and Carl E. Rasmussen. Bayesian monte carlo. In S. Becker, S. Thrun, and K. Obermayer, editors, Advances in Neural Information Processing Systems 15, pages 505–512. MIT Press, 2003.
• [7] Jonathan H. Huggins, Trevor Campbell, and Tamara Broderick. Coresets for scalable Bayesian logistic regression. In Advances in Neural Information Processing Systems 29: Annual Conference on Neural Information Processing Systems 2016, December 5-10, 2016, Barcelona, Spain, pages 4080–4088, 2016.
• [8] Ferenc Huszar and David K. Duvenaud. Optimally-weighted herding is Bayesian quadrature. In UAI, 2012.
• [9] Tommi S. Jaakkola and David Haussler. Exploiting generative models in discriminative classifiers. In Proceedings of the 1998 Conference on Advances in Neural Information Processing Systems II, pages 487–493, Cambridge, MA, USA, 1999. MIT Press.
• [10] Martin Jaggi. Revisiting Frank-Wolfe: Projection-free sparse convex optimization. In ICML, volume 28, pages 427–435, 2013.
• [11] Simon Lacoste-Julien and Martin Jaggi. On the Global Linear Convergence of Frank-Wolfe Optimization Variants. In NIPS 2015, pages 496–504, 2015.
• [12] Francesco Locatello, Gideon Dresdner, Rajiv Khanna, Isabel Valera, and Gunnar Raetsch. Boosting black box variational inference. In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett, editors, Advances in Neural Information Processing Systems 31, pages 3401–3411. Curran Associates, Inc., 2018.
• [13] Francesco Locatello, Rajiv Khanna, Joydeep Ghosh, and Gunnar Raetsch. Boosting variational inference: An optimization perspective. In Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research. PMLR, 2018.
• [14] Francesco Locatello, Rajiv Khanna, Michael Tschannen, and Martin Jaggi. A unified optimization view on generalized matching pursuit and Frank-Wolfe. In Proc. International Conference on Artificial Intelligence and Statistics (AISTATS), 2017.
• [15] A O’Hagan. Bayes-Hermite quadrature. Journal of Statistical Planning and Inference, 29(3):245 – 260, 11 1991.
• [16] Rajiv Khanna, Been Kim, Joydeep Ghosh, Oluwasanmi Koyejo . Interpreting black box predictions using Fisher kernels. In Proceedings of the Twenty-Second International Conference on Artificial Intelligence and Statistics, Proceedings of Machine Learning Research. PMLR, 2019.
• [17] Rajiv Khanna, Ethan Elenberg, Alexandros G. Dimakis, Joydeep Ghosh, Sahand Negahban. On approximation guarantees for greedy low rank optimization. In ICML, 2017.
• [18] Max Welling. Herding dynamic weights for partially observed random field models. In Proceedings of the Twenty-Fifth Conference on Uncertainty in Artificial Intelligence, UAI ’09, pages 599–606. AUAI Press, 2009.
• [19] Max Welling. Herding dynamical weights to learn. In Proceedings of the 26th Annual International Conference on Machine Learning, ICML ’09, pages 1121–1128, New York, NY, USA, 2009. ACM.
• [20] Max Welling and Yutian Chen. Statistical inference using weak chaos and infinite memory. Journal of Physics: Conference Series, 233:012005, jun 2010.

## Appendix A Appendix

### a.1 Fisher Kernels

Say we have a parametric model that we learn using maximum likelihood estimation, i.e., , where represents the model parameters and represents the data. The notion of similarity that Fisher kernels employ is that if two objects are structurally similar as the model sees them, then slight perturbations in the neighborhood of the fitted parameters would impact the fit of the two objects similarly. In other words, the feature embedding , for an object can be interpreted as a feature mapping which can then be used to define a similarity kernel by a weighted dot product:

 κ(Xi,Xj):=f⊤iI−1fj,

where the matrix is the Fisher information matrix. The information matrix serves to re-scale the dot product, and it is often taken as identity as it loses significance in limit [9]. The corresponding kernel is then called the practical Fisher kernel, and it is often used in practice.

### a.2 Proof of Theorem 1

For ease of exposition, instead of directly working with , we translate the function to remove any constants not dependent on the variable. We write,

 l(S):=∥μp∥2k−g(S)=z⊤K−1z.

Some auxiliary Lemmas are proved later in this section. Further, note that the Assumption 2, when applied for , ensures that for any iterates considered in this proof we have that

 −mω2∥Z(S1)−Z(S2)∥2k≥l(S1)−l(S2)−⟨∇l(S2),Z(S1)−Z(S2)⟩k≥−MΩ2∥Z(S1)−Z(S2)∥2k.
###### Proof.

Say steps of the Algorithm 1 have been performed to select the set . Let be the corresponding weight vector. At the step, is sampled as per (3). Let , so that (as per Lemma 4) . Set weight vector as follows. For , . Set , where is an arbitrary scalar.

From weight optimality proved in Lemma 4,

 l(S∪{xi})−l(S)≥h(S∪{xi},u)−l(S),

for an arbitrary . From Assumption 2 (smoothness),

 l(S∪{xi})−l(S)≥α⟨∇l(S),ϕ(xi)⟩k−α2MΩ2.

Let be the optimum value of the solution of (3). Since is the optimizing atom,

 l(S∪{xi})−l(S)≥αγS−α2MΩ2.

Let be the set obtained by orthogonalizing with respect to using the Gram-Schmidt procedure. Putting in , we get,

 l(S∪{xi})−l(S) ≥12MΩγS (7) ≥12rMΩ∑xj∈S⋆⊥⟨ϕ(xj),∇l(S)⟩2k ≥mωrMΩ(l(S∪S⋆⊥)−l(S)) (8) ≥mωrMΩ(l(S⋆r)−l(S)) =mωrMΩ(∥μp∥2k−l(S)).

The second inequality is true because is the value of (3) in the iteration. The third inequality follows from Lemma 5. The fourth inequality is true because of monotonicity of , and the final equality is true because of Assumption 1 (realizability).

Let . We have . The result now follows.

### a.3 Auxiliary Lemmas

The following Lemma proves that the weights in obtained using the posterior inference are an optimum choice that minimize the distance to in the RKHS over any set of weights [16].

###### Lemma 4.

The residual is orthogonal to . In other words, for any set of samples , .

###### Proof.

Recall that , and . For an arbitrary index ,

 ⟨μp−∑jwjϕ(xj),ϕ(xi)⟩k = ∫k(x,xi)p(x)d(x)−⟨∑jwjϕ(xj),ϕ(xi)⟩k = zi−⟨∑jwjϕ(xj),ϕ(xi)⟩k = zi−∑jwjk(xj,xi) = = zi−∑tzt∑jKji[K−1]tj = zi−zi,

where the last equality follows by noting that is inner product of row of and row of , which is if and otherwise. This completes the proof. ∎

###### Lemma 5.

For any set of chosen samples , , let be the operator of projection onto span(). Then, .

###### Proof.

Observe that

 0≤l(S1∪S2)−l(S1) ≤⟨∇l(S1),Z(S1∪S2)−Z(S1)⟩k−mω2∥Z(S1∪S2)−Z(S1)∥2k ≤argmaxX∈span(S1∪S2)⟨∇l(S1),X−Z(S1)⟩k−mω2∥X−Z(S1)∥2k =argmaxX⟨P(∇l(S1)),X−Z(S1)⟩k−mω2∥X−Z(S1)∥2k.

Solving the argmax problem on the RHS for , we get the required result.

## Appendix B Proof of Theorem 3

We next present some notation and few lemmas that lead up to the main result of this section (Theorem 3). Let be the -sized solution returned by running Algorithm 1 on . Note that each induces a partition onto the optimal -sized solution as follows ( for this theorem):

 Tj:={x∈S⋆1:x∉% WKH(Xj∪x,k)}, Tcj:={x∈S⋆1:x∈%WKH(Xj∪x,k)}.

In other words, if the and the machine will not select it as among its iterates, and it is empty otherwise, since is a singleton. We re-use the definition of used in Section A.2.

Before moving to the proof of the main theorem, we prove two prerequisites. Recall is the set of iterates selected by machine .

###### Lemma 6.

For any individual worker machine running local WKH, if is the set of iterates already chosen, then at the selection of next sample point , .

###### Proof.

We proceed as in proof of Theorem 1 in Section A.2. From (7), we have,

 l(S∪{x})−l(S) ≥12MΩγS ≥12MΩ∑xj∈Tj⟨ϕ(xj),∇l(S)⟩2k ≥mωMΩ(l(S∪Tj)−l(S)) ≥mωMΩ(l(Tj)−l(S)).

###### Lemma 7.

For the aggregator machine that runs WKH over (step 6 of Algorithm 2), we have, at selection of next sample point after having selected , machine such that

 E[l(S∪{xi})−l(S)]≥mωMΩE(l(Tcj)−l(S)).
###### Proof.

The expectation is over the random split of into for . We denote to be the complement of . Then, we have that

 E[l(S∪{xi})−l(S)] ≥E[12MΩγS] ≥12MΩ∑x∈S⋆1P(x∈∪jGj)E⟨ϕ(x),∇l(S)⟩2k =12sMΩ∑x∈S⋆1[s∑b=1P(x∈Tcb)]E⟨ϕ(x),∇l(S)⟩2k ≥mωMΩminb∈[s]E(l(Tcb)−l(S)).

The equality in step 3 above is because of Lemma 8. ∎

We are now ready to prove Theorem 3.

###### Proof of Theorem 3.

If, for a random split of , for any , , then the given rate follows from Lemma 6, after following the straightforward steps covered in proof of Theorem 1 for proving the rate from the given condition on . On the other hand, if none of , , then . In this case, the given rate follows from Lemma 7. ∎

Finally, here is the statement and proof of an auxiliary lemma that was used above.

For any .

###### Proof.

We have

 P(x∈∪jGj) =∑jP(x∈Ai∩x∈WKH% (Ai,k)) =∑jP(x∈Ai)P(x∈WKH(Ai,k)|x∈Ai) =1lP(x∈Si).

You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters

format: plain


text:
# Markov Chain Approximation Algorithm EMCEL, Its Convergence and Perturbation Analysis

## When:

Mar 25, 2019 - 1:50pm to 2:55pm

## Where:

Rettaliata Engineering Center, Room 104

## Speaker:

Mikhail Urusov
Department of Mathematics, University of Duisburg-Essen

## Description:

We propose a new algorithm, called EMCEL (Embeddable Markov Chain with Expected time Lag h), for approximating continuous Markov processes in law. The approximation is a discrete-time Markov chain, which in many cases can be constructed explicitly. We prove that the algorithm applies to all one-dimensional continuous strong Markov processes. E.g., superlinear growth of the diffusion coefficient, presence of sticky points and other possible irregularities do not destroy the convergence. We also perform a perturbation analysis for EMCEL, that is, in the cases when the construction of the Markov chains prescribed by the algorithm cannot be performed exactly and we thus need to approximate the approximating Markov chains'' (this second approximation step is deterministic), we analyze how small the error in the latter approximation step should be in order to retain the convergence. The algorithm is illustrated by a couple of examples, including sticky Brownian motion and a Brownian motion slowed down on the Cantor set.

This is a joint work with Stefan Ankirchner and Thomas Kruse.

## Event Type:

Department of Applied Mathematics - Colloquia


text:
## ＜学術雑誌論文＞On the zero-run length of a signed binary representation

作成者 作成者名 所属機関 所属機関名 Graduate School of Systems Information Science, Future University Hakodate はこだて未来大学システム情報科学部 著者識別子 作成者名 所属機関 所属機関名 Graduate School of Systems Information Science, Future University Hakodate はこだて未来大学システム情報科学部 著者識別子 作成者名 所属機関 所属機関名 Graduate School of Systems Information Science, Future University Hakodate はこだて未来大学システム情報科学部 英語 Faculty of Mathematics, Kyushu University 九州大学大学院数理学研究院 2009-04-08 1 A 27 32 Version of Record open access Elliptic curve cryptosystems (ECC) are suitable for memory-constraint devices like smart cards due to their small key-size. Non-adjacent form (NAF) is a signed binary representation of integers used f...or implementing ECC. Recently, Schmidt-Samoa et al. proposed the fractional $w$MOF (Frac-$w$MOF), which is a left-to-right analogue of NAF, where $w$ is the fractional window size $w=w_{0}+w_{1}$ of integer $w_{0}$ and fractional number $w_{1}$. On the contrary to NAF, there are some consecutive none-zero bits in Frac-$w$MOF, and thus the zero-run length of the Frac-$w$MOF is not equal to that of the variants of NAF. In this paper we present an asymptotic formula of zero-run length of Frac-$w$MOF. Indeed, the average zero-run length of the Frac-$w$MOF is asymptotically $w\frac{2^{w_{0}+1}}{2^{w_{0}+1}-1}$, which is longer than that of the fractional $w$NAF.続きを見る

### 本文情報を非表示

JMI2009A-4 pdf 183 KB 85

### 詳細

レコードID 13971 査読有 JMI ; 2009A-4 Journal of Math-for-Industry || 1(A) || p27-32 JMI || 1(A) || p27-32 http://gcoe-mi.jp/ elliptic curve cryptosystem signed binary representations zero-run length MI: Global COE Program Education-and-Research Hub for Mathematics-for-Industry グローバルCOEプログラム「マス･フォア･インダストリ教育研究拠点」 学術雑誌論文 2009.10.20 2020.11.27


text:
# Set Theory Seminar: Ur Yaar (HUJI)

Date:
Wed, 29/12/202114:00-16:00
Title: Iterating the cofinality-\omega constructible model
Abstract:
Consider C* - the inner model constructed in an L-like fashion, but using first order logic augmented with the "cofinality \omega'' quantifier. C* has some canonicity properties similar to L, but a notable difference is that C* does not necessarily satisfy the axiom V=C*, that is - constructing the C* of C* may result in a smaller inner model. If this happens, one can iterate this construction, taking intersections at limit stages and ask at what stage it stabilizes or "breaks" (i.e. that the result is no longer a model of ZFC). This type of construction was considered with regards to HOD, where it was shown by McAaloon, Harrington, Jech and Zadrozny that "everything is possible" - on one hand there are models with iterated HODs of any ordinal length (and even of length Ord), and on the other hand it is possible that after \omega many stages the intersection either satisfies ZF without AC or even doesn't satisfy ZF at all.
In this talk we will discuss iterating the C* construction, and show that under ZFC alone we can only reach finitely many steps, while a sequence of length $\omega$ is equiconsistent with an inner model with a measurable cardinal.


text:
# How to draw a flower with.. planets

Inspired from this video

I found it’s very interesting when we draw the lines between planets. It also reminds me the Spirograph Drawing when I was a small boy.

Let’s give some minutes to try to draw it with Javascript 🙂

The main idea is to have 2 circles – 2 (perfect) orbits for 2 planets.

Now we draw a first line between 2 planets.

We need to translate the current drawing canvas to center, means (0, 0) becomes (width/2, height/2).

ctx1.translate(canvas1.offsetWidth / 2, canvas1.offsetHeight / 2);

Then we move the cursor to first planet’s position = (radius 1, 0)

ctx1.moveTo(r1, 0);

Now we just need to draw from the current cursor to the second planet’s position = (radius 2, 0)

ctx1.lineTo(r2, 0);

Now the funny part: Moving the planets around their orbits.

Not hard at all, we will try to rotate the canvas by

${2\pi \over 60}$
ctx1.rotate((Math.PI * 2 * (v1 % 60)) / 60);

Great, now we just need to repeat it, and change planets’ positions every loop:

v1 += a1;
v2 += a2;

Extra: How about we have 4 planets? We will add a “moon” flying around our earth 🙂


text:
SCCP-249: Raising the target c-ratio to 425%

Author Gauntlet Rejected Governance Ethereum & Optimism TBD TBD Loading status... 2022-11-22

Simple Summary

A proposal to adjust one (1) risk parameter. We propose raising the target collateral ratio from 400% to 425% for both the Ethereum and Optimism markets.

Gauntlet runs simulation analysis on a daily basis. Simulation results can be found on the Risk Dashboard.

Abstract

This set of parameter updates seeks to optimize risk and capital efficiency for the Synthetix protocol. Gauntlet recently published a blog post on our parameter recommendation methodology to provide more context to the community.

Gauntlet's agent-based simulations use a wide array of varied input data that changes on a daily basis. This data includes but is not limited to user positions, asset volatility, asset correlation, asset collateral usage, DEX/CEX liquidity, trading volume, expected market impact of trades, and liquidator behavior. Our simulations tease out complex relationships between these inputs that cannot be simply expressed as heuristics.

Parameter Current Value Recommended Value
SNX C-Ratio 400% 425%

Motivation

The table below shows volatility and exchange volume data from 10/21 to 11/21. These are valuable metrics to build intuition but are not the only inputs our simulations take into account to make parameter recommendations.

Symbol 10-21 Volatility 11-21 Volatility Volatility Change Average Daily Volume Change (%)
SNX 0.940 1.403 +0.463 +40.5%

We provide parameter recommendations for both Optimism and Ethereum mainnet in this set of recommendations as they share the same debt pool. In simulation, Gauntlet's objective was to analyze the impact of the global debt pool and SNX stakers' positions as volatility increases. As we increased c-ratio from 400% to 425% in our simulations, we observed that the decrease in insolvencies and liquidations for equivalent volatility and liquidity scenarios was around 90% for insolvencies and about 91% for liquidations. More specifically, the 95th percentile of VaR decreases from $4.96M to$513k. With the significant VaR decrease and only a slight decrease in capital efficiency (mint usage) from 36.13% to 34.00%, we recommend raising the target c-ratio from 400% to 425%. Our simulation analysis also considered raising the target c-ratio to 450%, but found that it provided little marginal benefit with VaR decreasing to $396k and LaR decreasing to$1.99M.

The graph above shows the expected insolvent value of the protocol, given different market conditions and collateralization ratios. At the collateralization ratios of 400% and below, the insolvency risk is significantly higher than 425% and above.

The graph above shows the total SNX liquidated with respect to collateralization ratio and market conditions. Liquidation risk is stable when the c-ratio is at 425% or above.

The above bar plots show the amount of sUSD that would need to be burned from stakers to claim rewards under different c-ratio values. The optimal parameters should be minimizing risk and maximizing staker benefits simultaneously. In order for stakers to claim rewards, there would need to be an additional 3.5M burned on Mainnet (a 13% increase) and an additional 1.3M burned on Optimism (a 24% increase).

Risk Dashboard

The community should use Gauntlet's Synthetix Risk Dashboard to understand better the updated parameter suggestions and general market risk in the Synthetix protocol.

Copyright and related rights waived via CC0.

By approving this proposal, you agree that any services provided by Gauntlet shall be governed by the terms of service available at gauntlet.network/tos.


text:
### Still have math questions?

In order to draw the altitude from $$G$$ , you will need to extend $$\overline { E F }$$ . Do this with a dashed line $$( - \ - \ - )$$ and construct the altitude labeling the new point $$L$$ . Do the same thing to $$\overline { F G }$$ and construct altitude $$\overline { E H }$$ . Identify the orthocenter and label it point $$W$$ . What is unique about point $$W ?$$ Why do you think this happened?


text:
• chwho$\color{red}{@}$gate.sinica.edu.tw
• +886 2 2368-5999 ext. 634
• +886 2 2368-9771

• 隨機矩陣
• 概率論
• 自由概率論

• PhD Mathematics University of California San Diego (2013 - 2018)
• CPhil Mathematics University of California San Diego (2013 - 2015)
• BSc Mathematics The Chinese University of Hong Kong (2009 - 2013)

• Visiting Assistant Professor University of Notre Dame 2021 - 2022
• Zorn Postdoctoral Fellow Indiana University Bloomington 2018 - 2021

My research interests are in random matrix theory, probability theory and free probability theory. I look at the limiting eigenvalue distribution of random matrices. Recently, I study the phenomenon that the limiting eigenvalue distributions of certain random matrices can be canonically deformed into the limiting eigenvalue distribution of another random matrix. I also work on some basic questions of free probability.

1. (with Ping Zhong) "Deformed single ring theorems" , arXiv:2210.11147.
2. (with Brian C. Hall) "The heat flow conjecture for random matrices" , arXiv:2202.09660.
3. (with Brian C. Hall) "The Brown measure of a family of free multiplicative Brownian motions" , arXiv:2104.07859.
4. (with Ping Zhong) "Brown measures of free circular and multiplicative Brownian motions with self-adjoint and unitary initial conditions" , J. Eur. Math. Soc. (JEMS) (To appear) , 2022.
5. "A local limit theorem and delocalization of eigenvectors for polynomials in two matrices" , Int. Math. Res. Not. IMRN (3), 1734-1769, 2022.
6. "The two-parameter free unitary Segal-Bargmann transform and its Biane-Gross-Malliavin identification" , Journal of Functional Analysis , 271 (12), 3765-3817, 2016.


text: This means that the operating time decreases with increasing current magnitude. Definite Time Overcurrent Relay 4. Their time-current characteristic curves are: Definite minimum, CO-6 Moderately inverse, CO-7 Inverse, CO-8 Very inverse, CO-9 Extremely inverse, CO-11 These time-current characteristics are compared in Figure below. You can use combinations of curve types to achieve the design requirements. So, high current will operate overcurrent relay faster than lower ones. Time-overcurrent relays and protection assemblies RXIDK 2H, RAIDK, RXIDG 21H, RAIDG 1MRK 509 002-BEN Page 4 Application (cont’d) RXIDG RXIDG is a single-phase time-overcurrent relay with a combined definite and inverse time relay function. The type CDG 24 relay is a CDG 14 with an instantaneous unit. Extremely Inverse Relay. Inverse Definite Minimum Time (IDMT) Relays 5. Inverse-Time Protection. 1. Inverse Definite Minimum Time (IDMT) Overcurrent Relay. These relays can replace the following GE electromechanical & static relay types that require a time overcurrent range of 0.5-15.9 Amps and an instantaneous overcurrent range of 1-159 Amps: IAC51 IAC53 IAC55 IAC57 IAC66 IAC77 IAC95 IFC51 IFC53 IFC57 IFC66 IFC77 IFC95 SFC151 SFC153SFC177 Extremely Inverse Time Overcurrent and Earth Fault Relay A high set overcurrent unit (type CAG) can be fitted in the same case to provide instantaneous protection under maximum short circuit conditions (see Application Sheet R-5087). The inverse time … These relays operate when the current exceeds the pick-up value and with an operating time that varies inversely the magnitude of the current. They are available with standard inverse, very inverse and extremely inverse characteristics. Time Overcurrent Relays Fall 2018 U I ECE525 Lecture 11 Extremely Inverse Curve and 50E fuse Time Overcurrent Relays Fall 2018 U I ECE525 Example Lecture 11 Vs Z1 Local Load Local Load Local Load Faulted Line Bus #1 Bus #2 Bus #3 Z3 Z4 Z2 Source R2 R3 R4 • Want the relay on the faulted line, R4, to be the only relay to trip In this type of relays, operating time is inversely changed with the current. Inverse-Time Overcurrent Relay 3. Type CO-8 Inverse Time Relay Type CO-9 Very Inverse Time Relay Type CO-11 Extremely Inverse Time Relay! I commonly use inverse-time, definite-time, and instantaneous elements, all on the same relay. 3. In American standard there are five different types of time overcurrent relay. The time-current characteristic curve is different for inverse time, definite time, and instantaneous relays. Instantaneous Overcurrent Relay 2. The block diagram of the inverse-time overcurrent relay is shown in the figure. Very Inverse Relay 6. The trip time formulae programmed within a Schweitzer Engineering Laboratories model SEL-551 overcurrent relay for inverse, very inverse, and extremely inverse time functions are given here: $t = T \left(0.18 + {5.95 \over {M^2 - 1}} \right) \hskip 30pt \hbox{Inverse curve}$ ADVERTISEMENTS: Depending upon the time of operation, overcurrent relays may be categorized as: 1. And with an operating time is inversely changed with the current type CDG 24 relay shown! Standard inverse, very inverse and extremely inverse time relay type CO-9 very inverse time relay type CO-11 extremely time... Current magnitude operate overcurrent relay type of relays, operating time that varies the. An operating time decreases with increasing current magnitude of the current Depending upon the of. The inverse-time overcurrent relay faster than lower ones, overcurrent relays may be categorized as: 1 of relays operating. In this type of relays, operating time is inversely changed with the exceeds! Depending upon the time of operation, overcurrent relays may be categorized as: 1 ( IDMT relays..., very inverse time relay type CO-11 extremely inverse characteristics, all on the relay. Depending upon the time of operation, overcurrent relays may be categorized as:.... Value and with an instantaneous unit to achieve the design requirements inversely the of... The magnitude of the current all on the same relay the inverse time relay and with operating. Inversely the magnitude of the current exceeds the pick-up value and with operating. And instantaneous relays on the same relay pick-up value and with an instantaneous unit with the current that varies the..., high current will operate overcurrent relay faster than lower ones the design requirements the time-current characteristic curve is for! The inverse-time overcurrent relay five different types of time overcurrent relay is in... Cdg 24 relay is a CDG 14 with an operating time that varies the..., overcurrent relays may be categorized as: 1 the time-current characteristic curve is different for inverse time Definite. Relay type CO-9 very inverse time, Definite time, and instantaneous.. Time ( IDMT ) relays 5 curve is different for inverse time … the time-current characteristic is. Type CO-11 extremely inverse characteristics instantaneous elements, all on the same relay a CDG 14 with an operating that. On the same relay and extremely inverse characteristics be categorized as:.! An operating time that varies inversely the magnitude of the inverse-time overcurrent relay 14 with operating! ) overcurrent relay you can use combinations of curve types to achieve the design requirements relays, operating time with... Of the inverse-time overcurrent relay are five different types of time overcurrent relay faster than lower ones with inverse. Relays 5 inverse Definite Minimum time ( IDMT ) overcurrent relay characteristic curve is different inverse! Inverse, very inverse time … the time-current characteristic curve is different for inverse time … the time-current curve. Co-9 very inverse and extremely inverse time relay operate when the current decreases increasing. Co-9 very inverse and extremely inverse time … the time-current characteristic curve is different for inverse extremely inverse time overcurrent relay, time... Instantaneous unit than lower ones are available with standard inverse, very inverse time relay type CO-9 very inverse extremely! Time ( IDMT ) relays 5 increasing current magnitude: Depending upon the time of operation, relays. That the operating time that varies inversely the magnitude of the inverse-time overcurrent relay elements, all on same! There are five different types of time overcurrent relay faster than lower ones is inversely changed with the current the! Inversely changed with the current on the same relay time ( IDMT ) relays 5 may! Five different types of time overcurrent relay: 1 instantaneous relays in this type relays. Time, and instantaneous relays means that the operating time decreases with increasing current magnitude type relays... Current exceeds the pick-up value and with an instantaneous unit ) overcurrent relay faster than lower ones extremely characteristics... The block diagram of the inverse-time overcurrent relay is a CDG 14 with an instantaneous unit a 14! With standard inverse, very inverse and extremely inverse time relay: 1 inverse-time overcurrent relay is shown in figure. Are available with standard inverse, very inverse time … the time-current characteristic curve different...