Pointwise characterizations in generalized function algebras

P.F. Stadler: Landscapes and Their Correlation Functions
Most of the knowledge about landscapes has so far been derived using statistical methods, considering random models of landscapes rather than a single landscape. The distribution of local optima and the statistical characteristics of down-hill walks have been computed for the uncorrelated landscape of the random energy model 25, 26, 27], Furthermore, two one-parameter families of tunably rugged landscapes have been studied extensively: the Nk model and its variants 28, 29, 30, 16] and the p-spin models 31, 32, 33, 34]. Local optima of 2-spin models are considered in 35, 36, 37]. While the statistical approach is the natural one, e.g., in the physics of spin glasses, it seems to be rather contrived in evolutionary biology because it is by no means clear what a reasonable statistical model should look like, even if there is a computational procedure to model the landscapes of, say, RNA free energies. A theory of landscapes is based on three ingredients: we are given a nite, but very large set V of \con gurations" and a \ tness function" f : V ! IR. The third ingredient is a notion of neighborhood between the con gurations, which allows us to interpret V as the vertex set of a graph ?. We will refer to ? as the con guration space of the landscape f . Let us brie y discuss two examples here: Consider the set of RNA molecules of given chain length n. A particular molecule x can be represented as a string of length n taken from the alphabet fG, C, A, Ug; molecular biologists call this string the sequence of the RNA. The \ tness" function f is, for instance, the free energy of folding x into its secondary structure 15]. In silico the folding is done by an algorithm containing a large number of experimentally determined parameters 38]. In nature as well as in in vitro experiments variation is introduced by mutations, predominantly point-mutations. Neighboring sequences are thus those that di er in only a single position. The resulting graph is known as the sequence space 39, 40]. A very di erent example is the Travelling Salesman Problem. A salesman starts from his home city and visits exactly once each of the n cities on a given list, then he returns home. The con gurations are the possible tours, i.e., all permutations of the cities on the salesman's list. The numerical value assigned to a particular tour {2{

. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
. . . . . . . . . . . . . . . . .
另一方面梁的轴向应变的表达式也会因变形大小的不同而采用不同的表达式比如小变形时应变而当考虑大变形时可能采用的应变表达式就是进而得到的梁的振动方程将会是一个含有高度非线性项的偏微分方程组
非线性动力学入门
张新华
西安交通大学 工程力学系 2011 年 07 月
前 言
─1687 年，牛顿(Isaac Newton, 1643 ～ 1727)发表了《自然哲学之数学原 理》(Mathematical Principles of Natural Philosophy)，标志着经典力学(亦即牛 顿力学)的正式诞生。牛顿力学主要研究自由质点系的宏观运动规律。 ─1788 年，拉格朗日(Joseph Louis Lagrange, 1736 ～ 1813)发表了分析力 学教程(Analytical Mechanics)，标志着拉格朗日力学的诞生。Lagrange 力学属 于分析力学的主要内容之一，在位形空间中研究带有约束的质点系动力学。 ─1833 年，哈密尔顿(William Rowen Hamilton, 1805 ～ 1865)对 Lagrange 力学进行了改造，引进了相空间(2n 维空间)，对系统内在的对称性(辛对称， Symplectic)进行了刻画。狭义上的哈密尔顿力学只适用于保守系统，而广义 的哈密尔顿力学在适用于非保守系统。哈密尔顿力学也属于分析力学的主要 组成部分。在此后发展起来的量子力学中 Hamilton 力学发挥着巨大的作用。 目前在天体力学、计算 Hamilton 力学，量子力学，甚至弹性力学(即所谓的 辛弹性力学)中哈密尔顿力学依然发挥着重要作用。 ─1927 年，Birkhoff(George David Birkhoff, 1844 ～ 1944)发表了“动力系 统”(Dynamical Systems)，标志着 Birkhoff 动力学的正式问世。Birkhoff 动力 学建立了研究非完整力学的框架。 ─1892 ～ 1899， 彭加莱(Henri Poincaré, 1854 ～ 1912)发表了三卷本的“天 体力学中的新方法”(New Methods of Celestial Mechanics)，系统性地提出了 研究动力学系统的定性方法，即几何方法。经典力学的目标之一就是设法求 得系统的解析解，而 Poincaré意识到对于大多数非线性系统而言，求其解析 解是不可能的，而必须发展新的研究方法。他超越了他的时代，极富远见地 预测到了非线性系统混沌现象(系统的解对初始条件具有极端敏感依赖性)的 存在。更为重要的是，Poincaré开创了研究非线性动力系统的几何方法，当之 无愧地被誉为非线性科学之父，其影响是划时代的。 ─1892 年，李亚普诺夫(Aleksandr Mikhailovich Lyapunov, 1857 ～ 1918)在 他的博士论文“运动稳定性的一般问题”(General problem of the stability of motion )中，系统地探讨了非线性动力学系统的稳定性问题。他提出了两种研 究稳定性的方法：李亚普诺夫第一方法(间接方法)和李亚普诺夫第二方法(直 接方法)。他从代数角度出发，对动力学系统的研究开创了一个崭新的领域。 彭加莱与李亚普诺夫，前者从几何角度，后者从代数角度，开拓了非线 性科学的研究疆域和研究手段。 ─1963 年，Lorenz(Edward Norton Lorenz, 1917 ～ 2008)发表了“确定性 非周期流”(Deterministic Nonperiodic Flow)的论文，认为大气系统的性态对 初值极为敏感，从而导致准确的长期天气预报是不可能的。该文标志着人类 首次借助于计算机发现了混沌(Chaos)现象的存在。 ─1757 年，欧拉(Leonhard Euler, 1707 ～ 1783)发表了压杆稳定性的论 文，首次探讨了力学系统的分岔现象。作为分岔理论重要分支的突变理 论(Catastrophe Theory)则主要由法国数学家托姆(René Thom, 1923 ～ 2002)于 上个世纪 60 年代创立，由齐曼(Christopher Zeeman，1925 ～)在 70 年代大力 推广普及。 ─1834 年，英国的罗素(John Scott Russell, 1808 ～ 1882)骑着马在 Union 运河上散步时，发现了现在称之为孤立波(又称作孤波，Solitary wave)的 i

Electric Power Systems Research 84 (2012) 58–64Contents lists available at SciVerse ScienceDirectElectric Power SystemsResearchj o u r n a l h o m e p a g e :w w w.e l s e v i e r.c o m /l o c a t e /e p srDynamic equivalence by an optimal strategyJuan M.Ramirez a ,∗,Blanca V.Hernández a ,Rosa Elvira Correa b ,1a Centro de Investigación y de Estudios Avanzados del I.P.N.,Av del Bosque 1145.Col El Bajio,Zapopan,Jal.,45019,MexicobUniversidad Nacional de Colombia –Sede Medellín.Facultad de Minas.Carrera 80#65-223Bloque M8,114,Medellín,Colombiaa r t i c l ei n f oArticle history:Received 21May 2011Received in revised form 26September 2011Accepted 29September 2011Keywords:Power system dynamics Equivalent circuitPhasor measurement unit Stability analysisPower system stabilitya b s t r a c tDue to the curse of dimensionality,dynamic equivalence remains a computational tool that helps to analyze large amount of power systems’information.In this paper,a robust dynamic equivalence is proposed to reduce the computational burden and time consuming that the transient stability studies of large power systems represent.The technique is based on a multi-objective optimal formulation solved by a genetic algorithm.A simplification of the Mexican interconnected power system is tested.An index is used to assess the proximity between simulations carried out using the full and the reduced model.Likewise,it is assumed the use of information stemming from power measurements units (PMUs),which gives certainty to such information,and gives rise to better estimates.© 2011 Elsevier B.V. All rights reserved.1.IntroductionOne way to speed up the dynamic studies of currently inter-connected power systems without significant loss of accuracy is to reduce the size of the system model by means of dynamic equiva-lents.The dynamic equivalent is a simplified dynamic model used to replace an uninterested part,known as an external part,of a power system model.This replacement aims to reduce the dimension of the original model while the part of interest remains unchanged [1–6].The phrases “Internal system”(IS)and “external system”(ES)are used in this paper to describe the area in question,and the remaining regions,respectively.Boundary buses and tie lines can be defined in each IS or ES.It is usually intended to perform detailed studies in the IS.However,the ES is important to the extent where it affects IS analyses.The equivalent does not alter the transient behavior of the part of the system that is of concern and greatly reduces the dimen-sion of the network,reducing computational time and effort [4,7,8].The dynamic equivalent also can meet the accuracy in engineering,achieving effective,rapid and precise stability analysis and security controls for large-scale power system [4,8].However,the determi-∗Corresponding author.Tel.:+523337773600.E-mail addresses:jramirez@gdl.cinvestav.mx (J.M.Ramirez),bhernande@gdl.cinvestav.mx (B.V.Hernández),elvira.correa@ (R.E.Correa).1Tel.:+57314255140.nation of dynamic equivalents may also be a time consuming task,even if performed off-line.Moreover,several dynamic equivalents may be required to represent different operating conditions of the same system.Therefore,it is important to have computational tools that automate the procedure to evaluate the dynamic equivalent [7].Ordinarily,dynamic equivalents can be constructed follow-ing two distinct approaches:(i)reduction approach,and (ii)identification approach.The reduction approach is based on an elimination/aggregation of some components of the existing model [4,5,9].The two mostly found in the literature are known as modal reduction [6,10]and coherency based aggregation [2,11,12].The identification approach is based on either parametric or non-parametric identification [13,14].In this approach,the dynamic equivalent is determined from online measurements by adjusting an assumed model until its response matches with measurements.Concerning the capability of the model,the dynamic equivalent obtained from the reduction approach is considerably more reli-able and accurate than those set up by the identification approach,because it is determined from an exact model rather than an approximation based on measurements.However,the reduction-based equivalent requires a complete set of modeling data (e.g.model,parameters,and operating status)which is rarely avail-able in practice,in particular the generators’dynamic parameters [5,13,15,16].On the other hand,due to the lack of complete system data,and/or frequently variations of the parameters with time,the importance of estimation methods is revealed noticeably.Especially,on-line model correction aids for employing adaptive0378-7796/$–see front matter © 2011 Elsevier B.V. All rights reserved.doi:10.1016/j.epsr.2011.09.023J.M.Ramirez et al./Electric Power Systems Research84 (2012) 58–6459Fig.1.190-buses46-generators power system.controllers,power system stabilizers(PSS)or transient stability assessment.The capability of such methods has become serious rival of the old conventional methods(e.g.the coherency[11,12] and the modal[6,10]approaches).The equivalent estimation meth-ods have spread,because it can be estimated founded on data measured only on the boundary nodes between the study system and the external system.This way,without any need of informa-tion from the external system,estimation process tries to estimate a reduced order linear model,which is replaced for the external part.Evidently,estimation methods can be used,in presence of perfect data of the network as well to compute the equivalent by simulation and/or model order reduction[15].Sophisticated techniques have become interesting subject for researchers to solve identification problems since90s.For example, to obtain a dynamic equivalent of an external subsystem,an opti-mization problem has been solved by the Levenberg–Marquardt algorithm[17].Artificial neural networks(ANN)are the most prevalent method between these techniques because of its high inherent ability for modeling nonlinear systems,including power system dynamic equivalents[15,18–25].Power system real time control for security and stability has pro-moted the study of on-line dynamic equivalent technologies,which progresses in two directions.One is to improve the original off-line method.The mainstream approach is to obtain equivalent model based on typical operation modes and adjust equivalent parame-ters according to real time collected information[4].Distributed coordinate calculation based on real time information exchanging makes it possible to realize on-line equivalence of multi-area inter-connected power system in power market environment[4].Ourari et al.[26]developed the slow coherency theory based on the struc-ture preservation technique,and integrated dynamic equivalence into power system simulator Hypersim,verifying the feasibility of on-line computation from both computing time and accuracy[27].Prospects of phasor measurement technique based on global positioning system(GPS)applied in transient stability control of power system are introduced in Ref.[4].Using real data collected by phasor measurement unit(PMU),with the aid of GPS and high-speed communication network,online dynamic equivalent of interconnected power grid may be achieved[4].In this paper,the dynamic equivalence problem is formulated by two objective functions.An evolutionary optimization method based on genetic algorithms is used to solve the problem.2.PropositionThe main objective of this paper is the external system’s model order reduction of an electrical grid,preserving only the frontier nodes.That is,those nodes of the external system directly linked to nodes of the study system.At such frontier nodes,fictitious gen-erators are allocated.The external boundary is defined by the user. Basically,it is composed by a set of buses,which connect the exter-nal areas to the study system.There is not restriction about this set. Different operating conditions are taken into account.work reductionAfirst condition for an equivalence strategy is the steady state preservation on the reduced grid;this means basically a precise60J.M.Ramirez et al./Electric Power Systems Research 84 (2012) 58–64Fig.2.Proposed strategy’s flowchart.voltage calculation.In this paper,all nodes of the external system are eliminated,except the frontier nodes.By a load flow study,the complex power that should inject some fictitious generators at such nodes can be calculated.The nodal balance equation yields,j ∈Jp ij +Pg i +Pl i =0,∀i ∈I(1)where I is the set of frontier nodes;J is the set of nodes linked directly to the i th frontier node;p ij is the active power flowing from the i th to j th node;Pg i is the generation at the i th node;Pl i is the load at the i th node.Thus,the voltages for the reduced model become equal to those of the full one.For studies where unbalanced conditions are important,a similar procedure could be followed for the negative and zero equivalent sequences calculation.2.2.Studied systemThe power system shown in Fig.1depicts a reduced version of the Mexican interconnected power system.It encompasses 7regional systems,with a generation capacity of 54GW in 2004and an annual consumption level of 183.3TWh in 2005.The transmis-sion grid comprises a large 400/230kV system stretching from the southern border with Central America to its northern interconnec-tions with the US.The grids at the north and south of the country are long and sparsely interconnected transmission paths.The major load centers are concentrated on large metropolitan areas,mainly Mexico City in the central system,Guadalajara City in the western system,and Monterrey City in the northeastern system.The subsystem on the right of the dotted line is considered as the system under study.Thus,the subsystem on the left is the exter-nal one.There are five frontier nodes (86,140,142,148and 188)and six frontier lines (86–184,140–141,142–143,148–143(2)and 188–187).Thus,the equivalent electrical grid has five fictitious generators at nodes 86,140,142,148and 188.Transient stabil-ity models are employed for generators,equipped with a static excitation system;its formulation is described as follows,dıdt=ω−ω0(2)dωdt=1T j [Tm −Te −D (ω−ω0)](3)dE q dt=1T d 0 −E g −(x d −x d )i d +E fd(4)dE ddt =1T d 0−E d +(x q −xq )i q (5)dE fd dt=1T A−E fd +K A (V ref +V s −|V t |) (6)where ı(rad)and ω(rad/s)represent the rotor angular positionand angular velocity;E d (pu)and E q(pu)are the internal transient voltages of the synchronous generator;E fd (pu)is the excitationvoltage;i d (pu)and i q (pu)are the d -and q -axis currents;T d 0(s)and T q 0(s)are the d -and q -open-circuit transient time constants;x’d (pu)and x q(pu)are the d -and q -transient reactances;Tm (pu)and Te (pu)are the mechanical and electromagnetic nominal torque;Tj is the moment of inertia;D is the damping factor;K A and T A (s)are the system excitation gain and time constant;V ref is the voltage reference;V t is the terminal voltage;V s is the PSS’s output (if installed).The corresponding parameters are selected as typical [28].2.3.FormulationGiven some steady state operating point (#CASES )the following objective functions are defined,min f =[f 1f 2]f 1=#CASESop =1w 1opNg intk =1ωk ori (t )−ωkequiv (x,t )2(7)f 2=#CASESop =1w 2opNg intk =1Pe kori (t )−Pe k equiv (x,t )2(8)subject to:0=Ngen eqj =1H j −ni =1i ∈LH i(9)where ωk ori is the time behavior of the angular velocity of those generators in the original system,that will be preserved (Ng int),J.M.Ramirez et al./Electric Power Systems Research84 (2012) 58–6461Fig.3.Fitness assignment of NSGA-II in the two-objective space.Fig.4.Case1:from top to bottom(i)angular position37(referred to slack);(ii)angular speed28;(iii)electrical torque41,after a three-phase fault at bus172.62J.M.Ramirez et al./Electric Power Systems Research 84 (2012) 58–64Fig.5.Case 3:from top to bottom (i)angular position 40(referred to slack);(ii)angular speed 34;(iii)electrical torque 39,after a three-phase fault at bus 144.after a disturbance within the internal area;ωk equiv is the time behavior of the angular velocity of those generators in the equiv-alent system,after the same disturbance within the internal area;Pe k ori is the time behavior of the electrical power of those gen-erators of the original system within the internal area;Pe k equiv is the time behavior of the electrical power of those generators in the equivalent system within the internal area;H k is the k th genera-tor’s inertia;L is the set of generators that belong to the external system;Ngen eq is the number of equivalent generators [16,25].The set of voltages S ={V i ,V j ,...,V k |complex voltages stemming from PMUs }has been included in the solution.The main challenge in a multi-objective optimization environ-ment is to minimize the distance of the generated solutions to the Pareto set and to maximize the diversity of the developed Pareto set.A good Pareto set may be obtained by appropriate guiding of the search process through careful design of reproduction oper-ators and fitness assignment strategies.To obtain diversification special care has to be taken in the selection process.Special care is also to be taken to prevent non-dominated solutions from being lost.Elitism addresses the problem of losing good solutions dur-ing the optimization process.In this paper,the NSGA-II algorithm (Non dominated Sorting Genetic Algorithm-II)[29–31]is used to solve the formulation.The algorithm NSGA-II has demonstrated to exhibit a well performance;it is reliable and easy to handle.It uses elitism and a crowded comparison operator that keeps diver-sity without specifying any additional parameters.Pragmatically,it is also an efficient algorithm that has shown better results to solve optimization problems with multi-objective functions in a series of benchmark problems [31,32].There are some other meth-ods that may be used.For instance,it is possible to use at least two population-based non-Pareto evolutionary algorithms (EA)and two Pareto-based EAs:the Vector Evaluated Genetic Algorithm (VEGA)[36],an EA incorporating weighted-sum aggregation [33],the Niched Pareto Genetic Algorithm [34,35],and the Nondom-inated Sorting Genetic Algorithm (NSGA)[29–31];all but VEGA use fitness sharing to maintain a population distributed along the Pareto-optimal front.3.ResultsIn this case,the decision variables,x ,are eight parametersper each equivalent generator:{x d ,x d ,x q ,x q ,T d 0,T q 0,H,D }.In this paper,for five equivalent generators,there are 40parameters to be estimated.Likewise,in this case,a random change in the load of all buses gives rise to the transient behavior.A normal distribution with zero mean is utilized to generate the increment (decrement)in all buses.The variation is limited to a maximum of 50%.The disturbance lasts for 0.12s and then it is eliminated;the studied time is 2.0s.To attain more precise equivalence for severe operating conditions,this improvement could require load variations greater than 50%.However,this bound was used in all cases.Fig.2depicts a flowchart of the followed strategy to calculate an optimal solution.In this paper,three operating points are taken into account:(i)Case 1,the nominal case [37];(ii)Case 2,an increment of 40%in load and generation;(iii)Case 3,a decrement of 30%in load and gen-eration.To account for each operating condition into the objective functions,the same weighted factors have been utilized (w i =1/3),Eqs.(7)–(8).Table 1summarizes the estimated parameters for five equiv-alent generators,according to the two objective functions.TheJ.M.Ramirez et al./Electric Power Systems Research84 (2012) 58–6463 Table1Parameters of the equivalent under a maximum of50%in load variation.G equiv1G equiv2G equiv3G equiv4G equiv5f1f2f1f2f1f2f1f2f1f2x d0.1080.104 2.140 2.100 1.920 1.9100.1590.1790.3240.362x d0.1130.1130.8190.8690.3960.3960.08760.132 1.900 1.950T d010.7010.7039.6039.6018.8018.7011.5011.5011.7011.70x q0.7100.7500.5300.5190.2880.308 2.470 2.4600.7050.803x q0.3950.3850.8740.8050.9010.9260.9220.9530.8820.884T q0 4.950 4.82033.3033.40 4.280 4.10016.9016.9012.8012.90H 4.610 4.61022.2722.2747.2647.2669.2369.2333.9333.93D18.0318.40535.6536.6239.2239.141.9042.00705.1705.2 electromechanical modes associated to generators of the internalsystem are closely preserved.These generators arefictitious andbasically are useful to preserve some of the main interarea modesbetween the internal area and the external one[16].Thus,in order to avoid the identification of the equivalent gener-ators’parameters based on a specific disturbance,in this paper theuse of random changes in all the load buses is used.This will giverise to parameters valid for different fault locations.The allowedchange in the load(in this paper,50%)will result in a slight varia-tion of transient reactances.Further studies are required to assesssensitivities.Fig.3shows a typical Pareto front for this application.The NSGA-II runs on a Matlab platform and the convergence lastsfor3.25h for a population of200individuals and20generations.It is assumed that phasor measurement units(PMUs)areinstalled at specific buses(188,140,142,148,and86in the externalsystem,and141,143,145,and182in the internal system),whichbasically correspond to the frontier nodes.Likewise,it is assumedthat the precise voltages are known in these buses every time.Bythe inclusion of the PMUs,there is a noticeable improvement in thevoltages’information at the buses near them,due to the fact that itis assumed the PMUs’high precision.In this paper,the simulation results obtained by the full andthe reduced system are compared by a closeness measure,themean squared error(MSE).The goal of a signalfidelity measureis to compare two signals by providing a quantitative score thatdescribes the degree of similarity/fidelity or,conversely,the levelof error/distortion between ually,it is assumed that oneof the signals is a pristine original,while the other is distorted orcontaminated by errors[38].Suppose that z={z i|i=1,2,...,N}and y={y i|i=1,2,...,N}aretwofinite-length,discrete signals,where N is the number of signalsamples and z i and y i are the values of the i th samples in z and y,respectively.The MSE is defined by,MSE(z,y)=1NNi=1(z i−y i)2(10)Figs.4–5illustrate the transient behavior of some representa-tive signals after a three-phase fault at buses172(Fig.4)for theCase1;and144(Fig.5)for the Case3.Bus39is selected as theslack bus.Values in Table2show the corresponding MSE valuesfor the twelve generators of the internal system for each operat-ing case.Such values indicate a close relationship between the fulland reduced signals’behavior.In order to improve the equivalenceof a specific operating point,it is possible to weight it differentlythrough the factors w1–3,Eqs.(7)–(8).Table2shows the MSE’s val-ues when Case2has a higher weighting than Case1and Case3(w2=2/3,w1=w3=1/6).These values indicate that a closer agree-ment is attained between signals with the full and the reducedmodel.It is emphasized that the equivalence’s improvement couldrequire load variations greater than50%for off-nominal operatingconditions.4.ConclusionsUndoubtedly,the power system equivalents’calculationremains a useful strategy to handle the large amount of data,cal-culations,information and time,which represent the transientstability studies of modern power grids.The proposed approachis founded on a multi-objective formulation,solved by a geneticalgorithm,where the objective functions weight independentlyeach operating condition taken into account.The use of informationstemming from PMUs helps to improve the estimated equivalentgenerators’parameters.Results indicate that the strategy is ableto closely preserve the oscillating modes associated to the inter-nal system’s generators,under different operating conditions.Thatis due to the preservation of the machines’inertia.The use ofan index to measure the proximity between the signal’s behav-ior after a three-phase fault,indicates that good agreement isattained.In this paper,the same weighting factors have been usedto assess different operating conditions into the objective func-tions.However,depending on requirements,these factors can bemodified.The equivalence based on an optimal formulation assuresTable2MSE for Case2(Three-Phase Fault at Bus168).Angular position Angular speed Electrical powerw k=1/3w2=2/3,w1=w3=1/6w k=1/3w2=2/3,w1=w3=1/6w k=1/3w2=2/3,w1=w3=1/6 Gen39 2.13E−01 1.54E−01 6.36E−05 6.90E−059.40E−02 6.91E−02Gen28 1.79E−019.39E−02 3.43E−05 2.75E−05 4.09E−05 2.50E−05Gen29 1.12E−01 4.96E−02 4.84E−05 2.29E−059.25E−04 5.38E−04Gen30 1.01E−01 4.27E−02 3.31E−05 1.62E−059.75E−04 4.52E−04Gen32 2.80E−01 1.02E−01 4.34E−05 2.69E−05 4.26E−04 5.16E−04Gen33 2.97E−01 1.09E−01 4.81E−05 3.29E−059.85E−04 1.15E−04Gen34 1.70E−01 1.04E−01 4.03E−05 4.09E−059.19E−03 6.17E−03Gen37 1.68E−01 1.08E−01 5.23E−05 5.39E−05 1.12E−028.48E−03Gen38 1.53E−019.54E−02 4.47E−05 4.60E−05 1.29E−02 1.01E−02Gen40 1.91E−01 1.15E−01 5.27E−05 5.14E−05 2.26E−03 1.06E−03Gen41 1.93E−01 1.18E−01 5.52E−05 5.39E−05 2.50E−04 1.16E−04Gen42 1.68E−01 1.07E−01 4.17E−05 4.28E−05 5.23E−05 3.02E−0564J.M.Ramirez et al./Electric Power Systems Research84 (2012) 58–64proximity between the full and the reduced models.Closer prox-imity is reached if more stringent convergence’s parameters are defined,as well as additional objective functions,as line’s power flows,are included.References[1]P.Nagendra,S.H.nee Dey,S.Paul,An innovative technique to evaluate networkequivalent for voltage stability assessment in a widespread sub-grid system, Electr.Power Energy Syst.(33)(2011)737–744.[2]A.M.Miah,Study of a coherency-based simple dynamic equivalent for transientstability assessment,IET Gener.Transm.Distrib.5(4)(2011)405–416.[3]E.J.S.Pires de Souza,Stable equivalent models with minimum phase transferfunctions in the dynamic aggregation of voltage regulators,Electr.Power Syst.Res.81(2011)599–607.[4]Z.Y.Duan Yao,Z.Buhan,L.Junfang,W.Kai,Study of coherency-based dynamicequivalent method for power system with HVDC,in:Proceedings of the2010 International Conference on Electrical and Control Engineering,2010.[5]T.Singhavilai,O.Anaya-Lara,K.L.Lo,Identification of the dynamic equiva-lent of a power system,in:44th International Universities’Power Engineering Conference,2009.[6]J.M.Undrill,A.E.Turner,Construction of power system electromechanicalequivalents by modal analysis,IEEE Trans.PAS90(1971)2049–2059.[7]A.B.Almeida,R.R.Rui,J.G.C.da Silva,A software tool for the determinationof dynamic equivalents of power systems,in:Symposium-Bulk Power System Dynamics and Control–VIII(IREP),2010.[8]A.Akhavein,M.F.Firuzabad,R.Billinton,D.Farokhzad,Review of reductiontechniques in the determination of composite system adequacy equivalents, Electr.Power Syst.Res.80(2010)1385–1393.[9]U.D.Annakkage,N.C.Nair,A.M.Gole,V.Dinavahi,T.Noda,G.Hassan,A.Monti,Dynamic system equivalents:a survey of available techniques,in:IEEE PES General Meeting,2009.[10]B.Marinescu,B.Mallem,L.Rouco,Large-scale power system dynamic equiva-lents based on standard and border synchrony,IEEE Trans.Power Syst.25(4) (2010)1873–1882.[11]R.Podmore,A.Germond,Development of dynamic equivalents for transientstability studies,Final Report on EPRI Project RP763(1977).[12]E.J.S.Pires de Souza,Identification of coherent generators considering the elec-trical proximity for drastic dynamic equivalents,Electr.Power Syst.Res.78 (2008)1169–1174.[13]P.Ju,L.Q.Ni,F.Wu,Dynamic equivalents of power systems with onlinemeasurements.Part1:theory,IEE Proc.Gener.Transm.Distrib.151(2004) 175–178.[14]W.W.Price,D.N.Ewart,E.M.Gulachenski,R.F.Silva,Dynamic equivalents fromon-line measurements,IEEE Trans.PAS94(1975)1349–1357.[15]G.H.Shakouri,R.R.Hamid,Identification of a continuous time nonlinear statespace model for the external power system dynamic equivalent by neural net-works,Electr.Power Energy Syst.31(2009)334–344.[16]J.M.Ramirez,Obtaining dynamic equivalents through the minimizationof a lineflows function,Int.J.Electr.Power Energy Syst.21(1999) 365–373.[17]J.M.Ramirez,R.Garcia-Valle,An optimal power system model order reductiontechnique,Electr.Power Energy Syst.26(2004)493–500.[18]E.de Tuglie,L.Guida,F.Torelli,D.Lucarella,M.Pozzi,G.Vimercati,Identificationof dynamic voltage–current power system equivalents through artificial neural networks,Bulk power system dynamics and control–VI,2004.[19]A.M.Azmy,I.Erlich,Identification of dynamic equivalents for distributionpower networks using recurrent ANNs,in:IEEE Power System Conference and Exposition,2004,pp.348–353.[20]P.Sowa,A.M.Azmy,I.Erlich,Dynamic equivalents for calculation of powersystem restoration,in:APE’04Wisla,2004,pp.7–9.[21]A.H.M.A.Rahim,A.J.Al-Ramadhan,Dynamic equivalent of external power sys-tem and its parameter estimation through artificial neural networks,Electr.Power Energy Syst.24(2002)113–120.[22]A.M.Stankovic, A.T.Saric,osevic,Identification of nonparametricdynamic power system equivalents with artificial neural networks,IEEE Trans.Power Syst.18(4)(2003)1478–1486.[23]A.M.Stankovic,A.T.Saric,Transient power system analysis with measurement-based gray box and hybrid dynamic equivalents,IEEE Trans.Power Syst.19(1) (2004)455–462.[24]O.Yucra Lino,Robust recurrent neural network-based dynamic equivalencingin power system,in:IEEE Power System Conference and Exposition,2004,pp.1067–1068.[25]J.M.Ramirez,V.Benitez,Dynamic equivalents by RHONN,Electr.Power Com-pon.Syst.35(2007)377–391.[26]M.L.Ourari,L.A.Dessaint,V.Q.Do,Dynamic equivalent modeling of large powersystems using structure preservation technique,IEEE Trans.Power Syst.21(3) (2006)1284–1295.[27]M.L.Ourari,L.A.Dessaint,V.Q.Do,Integration of dynamic equivalents in hyper-sim power system simulator,in:IEEE PES General Meeting,2007,pp.1–6. [28]P.M.Anderson,A.A.Fouad,Power System Control and Stability.Appendix D,The Iowa State University Press,1977.[29]K.Deb,Evolutionary algorithms for multicriterion optimization in engineer-ing design,in:Evolutionary Algorithms in Engineering and Computer Science, 1999,pp.135–161(Chapter8).[30]N.Srinivas,K.Deb,Multiple objective optimization using nondominated sort-ing in genetic algorithms,put.2(3)(1994)221–248.[31]K.Deb,A.Pratap,S.Agarwal,T.Meyarivan,A fast and elitist multi-objectivegenetic algorithm:NSGA-II,IEEE put.6(2)(2002)182–197.[32]K.Deb,Multi-Objective Optimization using Evolutionary Algorithms,1st ed.,John Wiley&Sons(ASIA),Pte Ltd.,Singapore,2001.[33]P.Hajela,C.-Y.Lin,Genetic search strategies in multicriterion optimal design,Struct.Optim.4(1992)99–107.[34]Jeffrey Horn,N.Nafpliotis,Multiobjective optimization using the niched paretogenetic algorithm,IlliGAL Report93005,Illinois Genetic Algorithms Laboratory, University of Illinois,Urbana,Champaign,July,1993.[35]J.Horn,N.Nafpliotis,D.E.Goldberg,A niched pareto genetic algorithm formultiobjective optimization,in:Proceedings of the First IEEE Conference on Evolutionary Computation,IEEE World Congress on Computational Computa-tion,vol.1,82–87,Piscataway,NJ,IEEE Service Center,1994.[36]J.D.Schafer,Multiple objective optimization with vector evaluated geneticalgorithms,in:J.J.Grefenstette(Ed.),Proceedings of an International Confer-ence on Genetic Algorithms and their Applications,1985,pp.93–100.[37]A.R.Messina,J.M.Ramirez,J.M.Canedo,An investigation on the use of powersystem stabilizers for damping inter-area oscillations in longitudinal power systems,IEEE Trans.Power Syst.13(2)(1998)552–559.[38]W.Zhou,A.C.Bovik,Mean squared error:love it or leave it?A new look at signalfidelity measures,IEEE Signal Process.Mag.26(1)(2009)98–117.。

Kernels and Regularization on GraphsAlexander J.Smola1and Risi Kondor21Machine Learning Group,RSISEAustralian National UniversityCanberra,ACT0200,AustraliaAlex.Smola@.au2Department of Computer ScienceColumbia University1214Amsterdam Avenue,M.C.0401New York,NY10027,USArisi@Abstract.We introduce a family of kernels on graphs based on thenotion of regularization operators.This generalizes in a natural way thenotion of regularization and Greens functions,as commonly used forreal valued functions,to graphs.It turns out that diffusion kernels canbe found as a special case of our reasoning.We show that the class ofpositive,monotonically decreasing functions on the unit interval leads tokernels and corresponding regularization operators.1IntroductionThere has recently been a surge of interest in learning algorithms that operate on input spaces X other than R n,specifically,discrete input spaces,such as strings, graphs,trees,automata etc..Since kernel-based algorithms,such as Support Vector Machines,Gaussian Processes,Kernel PCA,etc.capture the structure of X via the kernel K:X×X→R,as long as we can define an appropriate kernel on our discrete input space,these algorithms can be imported wholesale, together with their error analysis,theoretical guarantees and empirical success.One of the most general representations of discrete metric spaces are graphs. Even if all we know about our input space are local pairwise similarities between points x i,x j∈X,distances(e.g shortest path length)on the graph induced by these similarities can give a useful,more global,sense of similarity between objects.In their work on Diffusion Kernels,Kondor and Lafferty[2002]gave a specific construction for a kernel capturing this structure.Belkin and Niyogi [2002]proposed an essentially equivalent construction in the context of approx-imating data lying on surfaces in a high dimensional embedding space,and in the context of leveraging information from unlabeled data.In this paper we put these earlier results into the more principled framework of Regularization Theory.We propose a family of regularization operators(equiv-alently,kernels)on graphs that include Diffusion Kernels as a special case,and show that this family encompasses all possible regularization operators invariant under permutations of the vertices in a particular sense.2Alexander Smola and Risi KondorOutline of the Paper:Section2introduces the concept of the graph Laplacian and relates it to the Laplace operator on real valued functions.Next we define an extended class of regularization operators and show why they have to be es-sentially a function of the Laplacian.An analogy to real valued Greens functions is established in Section3.3,and efficient methods for computing such functions are presented in Section4.We conclude with a discussion.2Laplace OperatorsAn undirected unweighted graph G consists of a set of vertices V numbered1to n,and a set of edges E(i.e.,pairs(i,j)where i,j∈V and(i,j)∈E⇔(j,i)∈E). We will sometimes write i∼j to denote that i and j are neighbors,i.e.(i,j)∈E. The adjacency matrix of G is an n×n real matrix W,with W ij=1if i∼j,and 0otherwise(by construction,W is symmetric and its diagonal entries are zero). These definitions and most of the following theory can trivially be extended toweighted graphs by allowing W ij∈[0,∞).Let D be an n×n diagonal matrix with D ii=jW ij.The Laplacian of Gis defined as L:=D−W and the Normalized Laplacian is˜L:=D−12LD−12= I−D−12W D−12.The following two theorems are well known results from spectral graph theory[Chung-Graham,1997]:Theorem1(Spectrum of˜L).˜L is a symmetric,positive semidefinite matrix, and its eigenvaluesλ1,λ2,...,λn satisfy0≤λi≤2.Furthermore,the number of eigenvalues equal to zero equals to the number of disjoint components in G.The bound on the spectrum follows directly from Gerschgorin’s Theorem.Theorem2(L and˜L for Regular Graphs).Now let G be a regular graph of degree d,that is,a graph in which every vertex has exactly d neighbors.ThenL=d I−W and˜L=I−1d W=1dL.Finally,W,L,˜L share the same eigenvectors{v i},where v i=λ−1iW v i=(d−λi)−1L v i=(1−d−1λi)−1˜L v i for all i.L and˜L can be regarded as linear operators on functions f:V→R,or,equiv-alently,on vectors f=(f1,f2,...,f n) .We could equally well have defined Lbyf,L f =f L f=−12i∼j(f i−f j)2for all f∈R n,(1)which readily generalizes to graphs with a countably infinite number of vertices.The Laplacian derives its name from its analogy with the familiar Laplacianoperator∆=∂2∂x21+∂2∂x22+...+∂2∂x2mon continuous spaces.Regarding(1)asinducing a semi-norm f L= f,L f on R n,the analogous expression for∆defined on a compact spaceΩisf ∆= f,∆f =Ωf(∆f)dω=Ω(∇f)·(∇f)dω.(2)Both(1)and(2)quantify how much f and f vary locally,or how“smooth”they are over their respective domains.Kernels and Regularization on Graphs3 More explicitly,whenΩ=R m,up to a constant,−L is exactly thefinite difference discretization of∆on a regular lattice:∆f(x)=mi=1∂2∂x2if≈mi=1∂∂x if(x+12e i)−∂∂x if(x−12e i)δ≈mi=1f(x+e i)+f(x−e i)−2f(x)δ2=1δ2mi=1(f x1,...,x i+1,...,x m+f x1,...,x i−1,...,x m−2f x1,...,x m)=−1δ2[L f]x1,...,x m,where e1,e2,...,e m is an orthogonal basis for R m normalized to e i =δ, the vertices of the lattice are at x=x1e1+...+x m e m with integer valuedcoordinates x i∈N,and f x1,x2,...,x m=f(x).Moreover,both the continuous and the dis-crete Laplacians are canonical operators on their respective domains,in the sense that they are invariant under certain natural transformations of the underlying space,and in this they are essentially unique.Regular grid in two dimensionsThe Laplace operator∆is the unique self-adjoint linear second order differ-ential operator invariant under transformations of the coordinate system under the action of the special orthogonal group SO m,i.e.invariant under rotations. This well known result can be seen by using Schur’s lemma and the fact that SO m is irreducible on R m.We now show a similar result for L.Here the permutation group plays a similar role to SO m.We need some additional definitions:denote by S n the group of permutations on{1,2,...,n}withπ∈S n being a specific permutation taking i∈{1,2,...n}toπ(i).The so-called defining representation of S n consists of n×n matricesΠπ,such that[Ππ]i,π(i)=1and all other entries ofΠπare zero. Theorem3(Permutation Invariant Linear Functions on Graphs).Let L be an n×n symmetric real matrix,linearly related to the n×n adjacency matrix W,i.e.L=T[W]for some linear operator L in a way invariant to permutations of vertices in the sense thatΠ πT[W]Ππ=TΠ πWΠπ(3)for anyπ∈S n.Then L is related to W by a linear combination of the follow-ing three operations:identity;row/column sums;overall sum;row/column sum restricted to the diagonal of L;overall sum restricted to the diagonal of W. Proof LetL i1i2=T[W]i1i2:=ni3=1ni4=1T i1i2i3i4W i3i4(4)with T∈R n4.Eq.(3)then implies Tπ(i1)π(i2)π(i3)π(i4)=T i1i2i3i4for anyπ∈S n.4Alexander Smola and Risi KondorThe indices of T can be partitioned by the equality relation on their values,e.g.(2,5,2,7)is of the partition type [13|2|4],since i 1=i 3,but i 2=i 1,i 4=i 1and i 2=i 4.The key observation is that under the action of the permutation group,elements of T with a given index partition structure are taken to elements with the same index partition structure,e.g.if i 1=i 3then π(i 1)=π(i 3)and if i 1=i 3,then π(i 1)=π(i 3).Furthermore,an element with a given index index partition structure can be mapped to any other element of T with the same index partition structure by a suitable choice of π.Hence,a necessary and sufficient condition for (4)is that all elements of T of a given index partition structure be equal.Therefore,T must be a linear combination of the following tensors (i.e.multilinear forms):A i 1i 2i 3i 4=1B [1,2]i 1i 2i 3i 4=δi 1i 2B [1,3]i 1i 2i 3i 4=δi 1i 3B [1,4]i 1i 2i 3i 4=δi 1i 4B [2,3]i 1i 2i 3i 4=δi 2i 3B [2,4]i 1i 2i 3i 4=δi 2i 4B [3,4]i 1i 2i 3i 4=δi 3i 4C [1,2,3]i 1i 2i 3i 4=δi 1i 2δi 2i 3C [2,3,4]i 1i 2i 3i 4=δi 2i 3δi 3i 4C [3,4,1]i 1i 2i 3i 4=δi 3i 4δi 4i 1C [4,1,2]i 1i 2i 3i 4=δi 4i 1δi 1i 2D [1,2][3,4]i 1i 2i 3i 4=δi 1i 2δi 3i 4D [1,3][2,4]i 1i 2i 3i 4=δi 1i 3δi 2i 4D [1,4][2,3]i 1i 2i 3i 4=δi 1i 4δi 2i 3E [1,2,3,4]i 1i 2i 3i 4=δi 1i 2δi 1i 3δi 1i 4.The tensor A puts the overall sum in each element of L ,while B [1,2]returns the the same restricted to the diagonal of L .Since W has vanishing diagonal,B [3,4],C [2,3,4],C [3,4,1],D [1,2][3,4]and E [1,2,3,4]produce zero.Without loss of generality we can therefore ignore them.By symmetry of W ,the pairs (B [1,3],B [1,4]),(B [2,3],B [2,4]),(C [1,2,3],C [4,1,2])have the same effect on W ,hence we can set the coefficient of the second member of each to zero.Furthermore,to enforce symmetry on L ,the coefficient of B [1,3]and B [2,3]must be the same (without loss of generality 1)and this will give the row/column sum matrix ( k W ik )+( k W kl ).Similarly,C [1,2,3]and C [4,1,2]must have the same coefficient and this will give the row/column sum restricted to the diagonal:δij [( k W ik )+( k W kl )].Finally,by symmetry of W ,D [1,3][2,4]and D [1,4][2,3]are both equivalent to the identity map.The various row/column sum and overall sum operations are uninteresting from a graph theory point of view,since they do not heed to the topology of the graph.Imposing the conditions that each row and column in L must sum to zero,we recover the graph Laplacian.Hence,up to a constant factor and trivial additive components,the graph Laplacian (or the normalized graph Laplacian if we wish to rescale by the number of edges per vertex)is the only “invariant”differential operator for given W (or its normalized counterpart ˜W ).Unless stated otherwise,all results below hold for both L and ˜L (albeit with a different spectrum)and we will,in the following,focus on ˜Ldue to the fact that its spectrum is contained in [0,2].Kernels and Regularization on Graphs5 3RegularizationThe fact that L induces a semi-norm on f which penalizes the changes between adjacent vertices,as described in(1),indicates that it may serve as a tool to design regularization operators.3.1Regularization via the Laplace OperatorWe begin with a brief overview of translation invariant regularization operators on continuous spaces and show how they can be interpreted as powers of∆.This will allow us to repeat the development almost verbatim with˜L(or L)instead.Some of the most successful regularization functionals on R n,leading to kernels such as the Gaussian RBF,can be written as[Smola et al.,1998]f,P f :=|˜f(ω)|2r( ω 2)dω= f,r(∆)f .(5)Here f∈L2(R n),˜f(ω)denotes the Fourier transform of f,r( ω 2)is a function penalizing frequency components|˜f(ω)|of f,typically increasing in ω 2,and finally,r(∆)is the extension of r to operators simply by applying r to the spectrum of∆[Dunford and Schwartz,1958]f,r(∆)f =if,ψi r(λi) ψi,fwhere{(ψi,λi)}is the eigensystem of∆.The last equality in(5)holds because applications of∆become multiplications by ω 2in Fourier space.Kernels are obtained by solving the self-consistency condition[Smola et al.,1998]k(x,·),P k(x ,·) =k(x,x ).(6) One can show that k(x,x )=κ(x−x ),whereκis equal to the inverse Fourier transform of r−1( ω 2).Several r functions have been known to yield good results.The two most popular are given below:r( ω 2)k(x,x )r(∆)Gaussian RBF expσ22ω 2exp−12σ2x−x 2∞i=0σ2ii!∆iLaplacian RBF1+σ2 ω 2exp−1σx−x1+σ2∆In summary,regularization according to(5)is carried out by penalizing˜f(ω) by a function of the Laplace operator.For many results in regularization theory one requires r( ω 2)→∞for ω 2→∞.3.2Regularization via the Graph LaplacianIn complete analogy to(5),we define a class of regularization functionals on graphs asf,P f := f,r(˜L)f .(7)6Alexander Smola and Risi KondorFig.1.Regularization function r (λ).From left to right:regularized Laplacian (σ2=1),diffusion process (σ2=1),one-step random walk (a =2),4-step random walk (a =2),inverse cosine.Here r (˜L )is understood as applying the scalar valued function r (λ)to the eigen-values of ˜L ,that is,r (˜L ):=m i =1r (λi )v i v i ,(8)where {(λi ,v i )}constitute the eigensystem of ˜L .The normalized graph Lapla-cian ˜Lis preferable to L ,since ˜L ’s spectrum is contained in [0,2].The obvious goal is to gain insight into what functions are appropriate choices for r .–From (1)we infer that v i with large λi correspond to rather uneven functions on the graph G .Consequently,they should be penalized more strongly than v i with small λi .Hence r (λ)should be monotonically increasing in λ.–Requiring that r (˜L) 0imposes the constraint r (λ)≥0for all λ∈[0,2].–Finally,we can limit ourselves to r (λ)expressible as power series,since the latter are dense in the space of C 0functions on bounded domains.In Section 3.5we will present additional motivation for the choice of r (λ)in the context of spectral graph theory and segmentation.As we shall see,the following functions are of particular interest:r (λ)=1+σ2λ(Regularized Laplacian)(9)r (λ)=exp σ2/2λ(Diffusion Process)(10)r (λ)=(aI −λ)−1with a ≥2(One-Step Random Walk)(11)r (λ)=(aI −λ)−p with a ≥2(p -Step Random Walk)(12)r (λ)=(cos λπ/4)−1(Inverse Cosine)(13)Figure 1shows the regularization behavior for the functions (9)-(13).3.3KernelsThe introduction of a regularization matrix P =r (˜L)allows us to define a Hilbert space H on R m via f,f H := f ,P f .We now show that H is a reproducing kernel Hilbert space.Kernels and Regularization on Graphs 7Theorem 4.Denote by P ∈R m ×m a (positive semidefinite)regularization ma-trix and denote by H the image of R m under P .Then H with dot product f,f H := f ,P f is a Reproducing Kernel Hilbert Space and its kernel is k (i,j )= P −1ij ,where P −1denotes the pseudo-inverse if P is not invertible.Proof Since P is a positive semidefinite matrix,we clearly have a Hilbert space on P R m .To show the reproducing property we need to prove thatf (i )= f,k (i,·) H .(14)Note that k (i,j )can take on at most m 2different values (since i,j ∈[1:m ]).In matrix notation (14)means that for all f ∈Hf (i )=f P K i,:for all i ⇐⇒f =f P K.(15)The latter holds if K =P −1and f ∈P R m ,which proves the claim.In other words,K is the Greens function of P ,just as in the continuous case.The notion of Greens functions on graphs was only recently introduced by Chung-Graham and Yau [2000]for L .The above theorem extended this idea to arbitrary regularization operators ˆr (˜L).Corollary 1.Denote by P =r (˜L )a regularization matrix,then the correspond-ing kernel is given by K =r −1(˜L ),where we take the pseudo-inverse wherever necessary.More specifically,if {(v i ,λi )}constitute the eigensystem of ˜L,we have K =mi =1r −1(λi )v i v i where we define 0−1≡0.(16)3.4Examples of KernelsBy virtue of Corollary 1we only need to take (9)-(13)and plug the definition of r (λ)into (16)to obtain formulae for computing K .This yields the following kernel matrices:K =(I +σ2˜L)−1(Regularized Laplacian)(17)K =exp(−σ2/2˜L)(Diffusion Process)(18)K =(aI −˜L)p with a ≥2(p -Step Random Walk)(19)K =cos ˜Lπ/4(Inverse Cosine)(20)Equation (18)corresponds to the diffusion kernel proposed by Kondor and Laf-ferty [2002],for which K (x,x )can be visualized as the quantity of some sub-stance that would accumulate at vertex x after a given amount of time if we injected the substance at vertex x and let it diffuse through the graph along the edges.Note that this involves matrix exponentiation defined via the limit K =exp(B )=lim n →∞(I +B/n )n as opposed to component-wise exponentiation K i,j =exp(B i,j ).8Alexander Smola and Risi KondorFig.2.Thefirst8eigenvectors of the normalized graph Laplacian corresponding to the graph drawn above.Each line attached to a vertex is proportional to the value of the corresponding eigenvector at the vertex.Positive values(red)point up and negative values(blue)point down.Note that the assignment of values becomes less and less uniform with increasing eigenvalue(i.e.from left to right).For(17)it is typically more efficient to deal with the inverse of K,as it avoids the costly inversion of the sparse matrix˜L.Such situations arise,e.g.,in Gaussian Process estimation,where K is the covariance matrix of a stochastic process[Williams,1999].Regarding(19),recall that(aI−˜L)p=((a−1)I+˜W)p is up to scaling terms equiv-alent to a p-step random walk on the graphwith random restarts(see Section A for de-tails).In this sense it is similar to the dif-fusion kernel.However,the fact that K in-volves only afinite number of products ofmatrices makes it much more attractive forpractical purposes.In particular,entries inK ij can be computed cheaply using the factthat˜L is a sparse matrix.A nearest neighbor graph.Finally,the inverse cosine kernel treats lower complexity functions almost equally,with a significant reduction in the upper end of the spectrum.Figure2 shows the leading eigenvectors of the graph drawn above and Figure3provide examples of some of the kernels discussed above.3.5Clustering and Spectral Graph TheoryWe could also have derived r(˜L)directly from spectral graph theory:the eigen-vectors of the graph Laplacian correspond to functions partitioning the graph into clusters,see e.g.,[Chung-Graham,1997,Shi and Malik,1997]and the ref-erences therein.In general,small eigenvalues have associated eigenvectors which vary little between adjacent vertices.Finding the smallest eigenvectors of˜L can be seen as a real-valued relaxation of the min-cut problem.3For instance,the smallest eigenvalue of˜L is0,its corresponding eigenvector is D121n with1n:=(1,...,1)∈R n.The second smallest eigenvalue/eigenvector pair,also often referred to as the Fiedler-vector,can be used to split the graph 3Only recently,algorithms based on the celebrated semidefinite relaxation of the min-cut problem by Goemans and Williamson[1995]have seen wider use[Torr,2003]in segmentation and clustering by use of spectral bundle methods.Kernels and Regularization on Graphs9Fig.3.Top:regularized graph Laplacian;Middle:diffusion kernel with σ=5,Bottom:4-step random walk kernel.Each figure displays K ij for fixed i .The value K ij at vertex i is denoted by a bold line.Note that only adjacent vertices to i bear significant value.into two distinct parts [Weiss,1999,Shi and Malik,1997],and further eigenvec-tors with larger eigenvalues have been used for more finely-grained partitions of the graph.See Figure 2for an example.Such a decomposition into functions of increasing complexity has very de-sirable properties:if we want to perform estimation on the graph,we will wish to bias the estimate towards functions which vary little over large homogeneous portions 4.Consequently,we have the following interpretation of f,f H .As-sume that f = i βi v i ,where {(v i ,λi )}is the eigensystem of ˜L.Then we can rewrite f,f H to yield f ,r (˜L )f = i βi v i , j r (λj )v j v j l βl v l = iβ2i r (λi ).(21)This means that the components of f which vary a lot over coherent clusters in the graph are penalized more strongly,whereas the portions of f ,which are essentially constant over clusters,are preferred.This is exactly what we want.3.6Approximate ComputationOften it is not necessary to know all values of the kernel (e.g.,if we only observe instances from a subset of all positions on the graph).There it would be wasteful to compute the full matrix r (L )−1explicitly,since such operations typically scale with O (n 3).Furthermore,for large n it is not desirable to compute K via (16),that is,by computing the eigensystem of ˜Land assembling K directly.4If we cannot assume a connection between the structure of the graph and the values of the function to be estimated on it,the entire concept of designing kernels on graphs obviously becomes meaningless.10Alexander Smola and Risi KondorInstead,we would like to take advantage of the fact that ˜L is sparse,and con-sequently any operation ˜Lαhas cost at most linear in the number of nonzero ele-ments of ˜L ,hence the cost is bounded by O (|E |+n ).Moreover,if d is the largest degree of the graph,then computing L p e i costs at most |E | p −1i =1(min(d +1,n ))ioperations:at each step the number of non-zeros in the rhs decreases by at most a factor of d +1.This means that as long as we can approximate K =r −1(˜L )by a low order polynomial,say ρ(˜L ):= N i =0βi ˜L i ,significant savings are possible.Note that we need not necessarily require a uniformly good approximation and put the main emphasis on the approximation for small λ.However,we need to ensure that ρ(˜L)is positive semidefinite.Diffusion Kernel:The fact that the series r −1(x )=exp(−βx )= ∞m =0(−β)m x m m !has alternating signs shows that the approximation error at r −1(x )is boundedby (2β)N +1(N +1)!,if we use N terms in the expansion (from Theorem 1we know that ˜L≤2).For instance,for β=1,10terms are sufficient to obtain an error of the order of 10−4.Variational Approximation:In general,if we want to approximate r −1(λ)on[0,2],we need to solve the L ∞([0,2])approximation problemminimize β, subject to N i =0βi λi −r −1(λ) ≤ ∀λ∈[0,2](22)Clearly,(22)is equivalent to minimizing sup ˜L ρ(˜L )−r−1(˜L ) ,since the matrix norm is determined by the largest eigenvalues,and we can find ˜Lsuch that the discrepancy between ρ(λ)and r −1(λ)is attained.Variational problems of this form have been studied in the literature,and their solution may provide much better approximations to r −1(λ)than a truncated power series expansion.4Products of GraphsAs we have already pointed out,it is very expensive to compute K for arbitrary ˆr and ˜L.For special types of graphs and regularization,however,significant computational savings can be made.4.1Factor GraphsThe work of this section is a direct extension of results by Ellis [2002]and Chung-Graham and Yau [2000],who study factor graphs to compute inverses of the graph Laplacian.Definition 1(Factor Graphs).Denote by (V,E )and (V ,E )the vertices V and edges E of two graphs,then the factor graph (V f ,E f ):=(V,E )⊗(V ,E )is defined as the graph where (i,i )∈V f if i ∈V and i ∈V ;and ((i,i ),(j,j ))∈E f if and only if either (i,j )∈E and i =j or (i ,j )∈E and i =j .Kernels and Regularization on Graphs 11For instance,the factor graph of two rings is a torus.The nice property of factor graphs is that we can compute the eigenvalues of the Laplacian on products very easily (see e.g.,Chung-Graham and Yau [2000]):Theorem 5(Eigenvalues of Factor Graphs).The eigenvalues and eigen-vectors of the normalized Laplacian for the factor graph between a regular graph of degree d with eigenvalues {λj }and a regular graph of degree d with eigenvalues {λ l }are of the form:λfact j,l =d d +d λj +d d +d λ l(23)and the eigenvectors satisfy e j,l(i,i )=e j i e l i ,where e j is an eigenvector of ˜L and e l is an eigenvector of ˜L.This allows us to apply Corollary 1to obtain an expansion of K asK =(r (L ))−1=j,l r −1(λjl )e j,l e j,l .(24)While providing an explicit recipe for the computation of K ij without the need to compute the full matrix K ,this still requires O (n 2)operations per entry,which may be more costly than what we want (here n is the number of vertices of the factor graph).Two methods for computing (24)become evident at this point:if r has a special structure,we may exploit this to decompose K into the products and sums of terms depending on one of the two graphs alone and pre-compute these expressions beforehand.Secondly,if one of the two terms in the expansion can be computed for a rather general class of values of r (x ),we can pre-compute this expansion and only carry out the remainder corresponding to (24)explicitly.4.2Product Decomposition of r (x )Central to our reasoning is the observation that for certain r (x ),the term 1r (a +b )can be expressed in terms of a product and sum of terms depending on a and b only.We assume that 1r (a +b )=M m =1ρn (a )˜ρn (b ).(25)In the following we will show that in such situations the kernels on factor graphs can be computed as an analogous combination of products and sums of kernel functions on the terms constituting the ingredients of the factor graph.Before we do so,we briefly check that many r (x )indeed satisfy this property.exp(−β(a +b ))=exp(−βa )exp(−βb )(26)(A −(a +b ))= A 2−a + A 2−b (27)(A −(a +b ))p =p n =0p n A 2−a n A 2−b p −n (28)cos (a +b )π4=cos aπ4cos bπ4−sin aπ4sin bπ4(29)12Alexander Smola and Risi KondorIn a nutshell,we will exploit the fact that for products of graphs the eigenvalues of the joint graph Laplacian can be written as the sum of the eigenvalues of the Laplacians of the constituent graphs.This way we can perform computations on ρn and˜ρn separately without the need to take the other part of the the product of graphs into account.Definek m(i,j):=l ρldλld+de l i e l j and˜k m(i ,j ):=l˜ρldλld+d˜e l i ˜e l j .(30)Then we have the following composition theorem:Theorem6.Denote by(V,E)and(V ,E )connected regular graphs of degrees d with m vertices(and d ,m respectively)and normalized graph Laplacians ˜L,˜L .Furthermore denote by r(x)a rational function with matrix-valued exten-sionˆr(X).In this case the kernel K corresponding to the regularization operator ˆr(L)on the product graph of(V,E)and(V ,E )is given byk((i,i ),(j,j ))=Mm=1k m(i,j)˜k m(i ,j )(31)Proof Plug the expansion of1r(a+b)as given by(25)into(24)and collect terms.From(26)we immediately obtain the corollary(see Kondor and Lafferty[2002]) that for diffusion processes on factor graphs the kernel on the factor graph is given by the product of kernels on the constituents,that is k((i,i ),(j,j ))= k(i,j)k (i ,j ).The kernels k m and˜k m can be computed either by using an analytic solution of the underlying factors of the graph or alternatively they can be computed numerically.If the total number of kernels k n is small in comparison to the number of possible coordinates this is still computationally beneficial.4.3Composition TheoremsIf no expansion as in(31)can be found,we may still be able to compute ker-nels by extending a reasoning from[Ellis,2002].More specifically,the following composition theorem allows us to accelerate the computation in many cases, whenever we can parameterize(ˆr(L+αI))−1in an efficient way.For this pur-pose we introduce two auxiliary functionsKα(i,j):=ˆrdd+dL+αdd+dI−1=lrdλl+αdd+d−1e l(i)e l(j)G α(i,j):=(L +αI)−1=l1λl+αe l(i)e l(j).(32)In some cases Kα(i,j)may be computed in closed form,thus obviating the need to perform expensive matrix inversion,e.g.,in the case where the underlying graph is a chain[Ellis,2002]and Kα=Gα.Kernels and Regularization on Graphs 13Theorem 7.Under the assumptions of Theorem 6we haveK ((j,j ),(l,l ))=12πi C K α(j,l )G −α(j ,l )dα= v K λv (j,l )e v j e v l (33)where C ⊂C is a contour of the C containing the poles of (V ,E )including 0.For practical purposes,the third term of (33)is more amenable to computation.Proof From (24)we haveK ((j,j ),(l,l ))= u,v r dλu +d λv d +d −1e u j e u l e v j e v l (34)=12πi C u r dλu +d αd +d −1e u j e u l v 1λv −αe v j e v l dαHere the second equalityfollows from the fact that the contour integral over a pole p yields C f (α)p −αdα=2πif (p ),and the claim is verified by checking thedefinitions of K αand G α.The last equality can be seen from (34)by splitting up the summation over u and v .5ConclusionsWe have shown that the canonical family of kernels on graphs are of the form of power series in the graph Laplacian.Equivalently,such kernels can be char-acterized by a real valued function of the eigenvalues of the Laplacian.Special cases include diffusion kernels,the regularized Laplacian kernel and p -step ran-dom walk kernels.We have developed the regularization theory of learning on graphs using such kernels and explored methods for efficiently computing and approximating the kernel matrix.Acknowledgments This work was supported by a grant of the ARC.The authors thank Eleazar Eskin,Patrick Haffner,Andrew Ng,Bob Williamson and S.V.N.Vishwanathan for helpful comments and suggestions.A Link AnalysisRather surprisingly,our approach to regularizing functions on graphs bears re-semblance to algorithms for scoring web pages such as PageRank [Page et al.,1998],HITS [Kleinberg,1999],and randomized HITS [Zheng et al.,2001].More specifically,the random walks on graphs used in all three algorithms and the stationary distributions arising from them are closely connected with the eigen-system of L and ˜Lrespectively.We begin with an analysis of PageRank.Given a set of web pages and links between them we construct a directed graph in such a way that pages correspond。