Efficient Analytical Derivatives of Rigid-Body Dynamics using Spatial Vector Algebra

An essential need for many model-based robot control algorithms is the ability to quickly and accurately compute partial derivatives of the equations of motion. State of the art approaches to this problem often use analytical methods based on the chain rule applied to existing dynamics algorithms. Although these methods are an improvement over finite differences in terms of accuracy, they are not always the most efficient. In this paper, we contribute new closed-form expressions for the first-order partial derivatives of inverse dynamics, leading to a recursive algorithm. The algorithm is benchmarked against chain-rule approaches in Fortran and against an existing algorithm from the Pinocchio library in C++. Tests consider computing the partial derivatives of inverse and forward dynamics for robots ranging from kinematic chains to humanoids and quadrupeds. Compared to the previous open-source Pinocchio implementation, our new analytical results uncover a key computational restructuring that enables efficiency gains. Speedups of up to 1.4x are reported for calculating the partial derivatives of inverse dynamics for the 50-dof Talos humanoid.

I. INTRODUCTION
Rigid-body dynamics models are widely used in the development of control algorithms for quadruped and humanoid robots, with recursive dynamics algorithms at the core of many real-time controllers. Example use cases include current approaches to controller design via optimization [1], [2]. Likewise, robotics libraries such as Crocoddyl [3] and Drake [4] require the partial derivatives of rigid-body dynamics with respect to the state and control variables. The calculation of these derivatives represents the majority of the CPU time required in many optimal control applications [5].
There are multiple general-purpose methods that can be applied to rigid-body dynamics for computing derivatives. Finite differences [7]- [9] are simple to implement and are trivially parallelized, but suffer in accuracy ( Table I). The complex-step method [12] is accurate to machine precision but requires all functions to be computed in the complex plane, leading to additional overhead. Cossette et al. [13] extended this method to matrix Lie groups, providing direct applicability to rigidbody dynamics. An approach that avoids the overhead of the complex step is automatic differentiation (AD), also called algorithmic differentiation [4], [10], [11]. AD calculates the partials of an algorithm by automating the accumulation of chain-rule expressions through operator overloading [10] or source code transformation [4], [11]. The first approach overloads basic data types to compute partials concurrently with all arithmetic. The second approach builds an expression graph from source code and augments the source directly to calculate the partial derivatives.
Without relying on AD, but similar to it, others have developed analytical methods that accumulate chain-rule expressions on top of existing algorithms. Examples include many existing descriptions of recursive algorithms for dynamics derivatives [5], [15], [16]. The derivation in Ref. [5] applies chain rule on the classical two-pass Recursive-Newton-Euler Algorithm (RNEA) for inverse dynamics (ID). The presented method has computational complexity O(N d) for both the forward and backward pass of the RNEA derivatives, where N is the number of bodies in the system and d is the depth of the kinematic connectivity tree. However, the algorithm in the C++ Pinocchio (2.6.0) code [6] associated with [5] is distinctly different from the algorithm presented in [5]. The Pinocchio code includes a more efficient O(N ) forward pass, coupled with an O(N d) backward pass.
Derivative algorithms that directly carry out chain-rule on an existing algorithm (using AD or by hand) may not always be the most efficient. Other special-purpose methods have been proposed to fully exploit the structure of the rigid-body equations of motion. Such approaches often analytically differentiate closed-form equations of motion and then design an algorithm to compute the result. An example is Ref. [17] where partial derivatives are considered of the mass matrix M (q), the Coriolis matrix C(q,q), and the gravity vector g(q) with respect to the state variables q andq, leading to an O(N 3 ) algorithm.
Ref. [18] gives a recursive method for partial derivatives of inverse and forward dynamics (FD) for serial kinematic chains with single-DoF joints. Their approach includes an O(N 2 ) algorithm for the partials of FD, and an O(N d) algorithm for the partials of ID.
The main contribution of this paper is to extend previous analytical results on partial derivatives of ID [18] to generic multi-DoF joint robots with a fixed or  1 Algo-2,3 of [5] provides a chain-rule method for ID partial derivatives. 2 Pinocchio code for ID partial derivatives accompanying [5] but different from chain-rule. 3 This paper provides a Special-Purpose Analytical Method that outperforms the state-of-the-art in terms of speed. floating base. We first derive closed-form expressions for the first-order partial derivatives of ID. To the best of authors' knowledge, these results represent the first of their kind for general rigid-body systems (fixed or floating base) with multi-DoF joints. The new expressions lead immediately to an algorithm of order O(N d) that naturally generalizes the one shown in [18]. The method is fundamentally different than the straight chain-rule approach presented in [5]. However, the distinct algorithm in the Pinocchio code [6] ends up being directly related to the computations required in our algorithm. Despite the similarity, our new closed-form expressions uncover a key restructuring that accelerates computations. We use the relationship between FD and ID [5], [18] in an efficient manner to ultimately calculate the partial derivatives of FD with complexity O(N 2 ).

II. DERIVATIVES OF RIGID-BODY DYNAMICS
Rigid-Body Dynamics: For a rigid-body system, the state variables are the configuration q and the generalized velocity vectorq, while the control variable is the generalized torque vector τ . The equation of motion, also called the Inverse Dynamics (ID), is given by where M ∈ R n×n is the mass matrix, C ∈ R n×n is a Coriolis matrix, and g ∈ R n is the vector of generalized gravitational forces, and n is the DoF of the system. For a fixed state, ID calculates τ for a givenq (Eq. 1), and Forward Dynamics (FD) computesq for a given τ : Toward supporting control applications (e.g., [1], [9]), the objective of this work is to compute the partial derivatives of FD with respect to the state (q,q) and control variables (τ ). The most efficient algorithms for calculating ID and FD are the O(N ) Recursive-Newton-Euler Algorithm (RNEA) and the Articulated-Body Algorithm (ABA) respectively [19].
Derivatives of Inverse and Forward Dynamics: As presented in [5], a direct approach to compute the partial derivatives of ID is to manually differentiate the RNEA algorithm via chain-rule, denoted as RNEACR (Table II). A similar chain rule approach for the derivatives of ABA is denoted as ABACR. Because the ABA is more computationally intensive than RNEA, a more efficient approach for derivatives of FD first applies RNEACR, and then uses the relationship [5], [18] where the variable u ∈ {q,q}. From Eq. 3, ∂F D/∂τ can be directly calculated as M −1 , enabling re-use in Eq. 5. This method for the partials of FD via RNEA Chain Rule is abbreviated as FDCR.
The state-of-the-art numerical method for computing partial derivatives of FD is provided in the Pinocchio package [6] (Pinocchio FD Derivs, Table II). The method first computes the derivatives of ID (Pinocchio ID Derivs, Table II). It then computes M −1 and uses Eq. 5 to obtain the derivatives of FD.
In the following sections, Spatial Vector Algebra (SVA) is reviewed for dynamics analysis, followed by a theoretical development of analytical expressions that lead to an algorithm for partial derivatives of ID. The new expressions result in algorithm with a key difference from the state-of-the-art [6], enabling a speedup over it.

III. SPATIAL VECTOR ALGEBRA (SVA)
Notation: Spatial vectors are 6D vectors that combine the linear and angular aspects of a rigid-body motion or net force [19]. Spatial vectors are denoted with lowercase bold letters (e.g., a), while matrices are denoted with capitalized bold letters (e.g., A). Motion vectors, such as velocity and acceleration, belong to a 6D vector space denoted M 6 . Spatial vectors are usually expressed in either the ground coordinate frame or a body coordinate (local coordinate) frame. For example, the spatial velocity k v k ∈ M 6 of a body k expressed in the body frame is given by is the angular velocity expressed in a coordinate frame  fixed to the body, while k v k ∈ R 3 is the linear velocity of the origin of the body frame. When expressed in the ground frame, the spatial velocity for body k is denoted as 0 v k . The vector is again composed of angular and linear velocity components. However, the linear velocity is associated with the body-fixed point on body k that is coincident with the origin of the ground frame 0. In this paper, when the frame used to express a spatial vector is omitted, the ground frame is assumed. Force-like vectors, such as force and momentum, belong to another 6D vector space F 6 , which is dual to M 6 . A spatial cross-product between two motion vectors (v,u), written as (v×)u, is given by (Eq. 6). This operation can be understood as providing the time rate of change of u, when u is moving with a spatial velocity v. For a Cartesian vector, the matrix ω× is the classical 3D cross product operator. A spatial cross-product between a motion and a force vector is written as (v× * )f , as defined in Eq. 7.
An operator × * is defined by swapping the order of the cross product, such that (f × * )v = (v× * )f [20]. Further introduction to SVA is provided in Appendix A.
Connectivity: An open-chain kinematic tree with serial or branched connectivity ( Fig. 1) is considered with N links connected by joints, each with up to 6 DoF. Body i's parent toward the root of the kinematic tree is denoted as λ(i). ν(i) denotes the set of bodies in the subtree rooted at body i, while ν(i) denotes the set of bodies in ν(i) excluding the body i . We define i j if body i is in the path from body j to the base. The spatial velocities of the neighbouring bodies in the tree are related by v i = v λ(i) + S iqi , where S i is the joint motion subspace matrix for joint i [19] anḋ q i the joint rates for joint i. For a single DoF revolute joint, S i ∈ R 6 , andq i is a scalar. For a 6-DoF freemotion joint, S i ∈ R 6×6 , whileq i ∈ R 6 . The velocity v i can also be written as the sum of joint velocities over predecessors as v i = l i S lql . The derivative of joint motion subspace matrix due to the axis changing with respect to local coordinates (commonly denotedS i [19]) is assumed to be zero. The quantityṠ i = v i × S i signifies the rate of change of S i due to the local coordinate system moving.

Dynamics:
The spatial equation of motion [19] is given for body k as: where f k is the net spatial force on body k, I k is its spatial inertia [19], and a k is its spatial acceleration. Instead of treating gravity as an external force, a common trick is to accelerate the base upwards opposite of the gravitational acceleration (a 0 = −a g ), providing the acceleration of body k as: For later use, we decompose a k into terms from joint accelerations, and terms from joint rates, according to In a similar fashion, the net spatial force on body k is then decomposed as: where η k = I k γ k (η ∈ F 6 ) is the spatial force on the body caused by joint accelerationsq, and ζ k = v k × * I k v k + I k ξ k (ζ∈ F 6 ) gives the Coriolis and centripetal forces on body k. From the formulation of the RNEA [19], Proceeding to consider the system overall, for any particular joint i, Eq. 1 can be written as follows.
Then, using τ i = S T i f C i , Eq. 9 and 11: where I C i denotes the composite rigid-body inertia of the sub-tree rooted at body i, given by IV. ANALYTICAL PARTIAL DERIVATIVES USING SVA In this section, closed-form expressions are derived for the derivatives of ID. The reader less interested in the derivation may skip to Eq. 30 and Eq. 32 as a summary.

A. Building Blocks
Several kinematic identities are given in Appendix C as a basis for the main derivation. We denote by n j the number of DoFs for joint j. The motion subspace matrix for the joint is then given as S j = s j,1 · · · s j,nj where each spatial vector s j,p gives the p-th free-mode of motion for joint j. With a slight liberty of notation, ∂/∂q j,p denotes an operator for the directional derivative along this p-th free mode of a joint. For example, considering the case with j i gives identity J1: which provides the rate of change in S i with respect to relative motion s j,p at joint j earlier in the chain. For revolute joints, these derivatives are just conventional derivatives with respect to joint angles. For multi-DoF joints such as a floating base, they are more formally Lie derivatives. Considering a configuration-dependent vector u, we denote by ∂u ∂q j = ∂u /∂q j,1 · · · ∂u /∂q j,n j the matrix of derivatives associated with joint j. To illustrate, we derive identity J7 for later use in the section.
Considering j i and using the definition of γ i : Using J1 and switching the order of the cross product: Collecting all DoFs of joint j, Eq. 80 becomes: Using the definition of γ l , and summing over l: The subsequent derivations employ these formulae and the identities from Appendix C to obtain the partials of the terms in Eq. 11. Since the derivations are quite lengthy, the presentation focuses on explaining the methodology. A full derivation is provided in Appendix D.
Partial Derivative of Eq. 12: We derive the partial derivative of [M (q)q] i with respect to q j for all i, j ∈ {1, . . . , N }. First, consider the case when j i. Using the product rule of differentiation in Eq. 12: Using identity J3 for the term , and cancelling terms: Upon summing over the index k, the final expression is: For the case i≺j, we have that: Using the identities J3, J4 and J7 and cancelling terms: Summing over the index k results in: where η C j collects joint-acceleration-dependent forces for the subtree, calculated as η C j = k j η k . For ease of implementation, we consider a single case where j i. The indices i and j are switched in Eq. 24, to get the expression for for the case j ≺ i. Therefore, Eq. 25 gives the two expressions formulated for the general case j i.
Partial Derivative of Eq. 13: Similarly, we use identities J2-J6 (Appendix C) to get the first-order partial derivatives of [C(q,q)q] i with respect to q j for the case j i as shown in Eq. 26, 27. In these equations, the composite of ζ C i for the subtree is defined as ζ C i = k i ζ k , and the matrix B k (v k , I k ) is a body-level Coriolis matrix [20]- [22] given by: Partial Derivative of the Gravity Term: Using the identities J3 and J4, for the case j i, the partial derivative of g i (Eq. 14) with respect to q j is: Summary for Partials of ID w.r.t. q: The partials of the individual components are now collected together. Terms in Eqs. 25, 26, 27, and 29 are added to get the total expressions for ∂τ ∂q for the case j i. The full derivation is provided in Appendix E.
The spatial quantitiesΨ j andΨ j are defined as:Ψ The partials of τ with respect toq depend only on the the Coriolis terms Cq. Using the identities J8 and J9 leads to expressions for the case when j i(see Appendix F for full derivation):

D. Algorithm for First-Order Partials of ID
Algorithm 1 returns the terms in Eq. 30 and Eq. 32 and has computational complexity O(N d). The method operates with all spatial vectors expressed in ground frame coordinates. The first pass is in the forward direction and goes from root to leaves calculating the spatial velocity, acceleration,Ṡ i ,Ψ i ,Ψ i , B i and the spatial force f i . When joint i has a single DoF,Ṡ i =Ψ i andS i =Ψ i and the steps are the same as in [18,Alg. 2].
The second pass progresses in the backward direction from leaves to the root. Index i takes all the values from 1 to N . The quantity ∂τ ∂q [ν(i), i] denotes the partial derivatives of τ for all bodies in the sub-tree of i with respect to q i , while ∂τ ∂q [i, ν(i)] is the partial derivative of τ for body i with respect to each of the bodies in the sub-tree of i, excluding body i. This backward pass again generalizes the one in [18,Alg. 2]. The differences are that the third term on line 15 can be dropped in the single-DoF case (since S T i (f C i × * )S i = 0 in that case), while lines 17-20 restructure a second inner loop in [18] into matrix multiplies.

A. Algorithm correctness
The partial derivatives of FD from ABACR, FDSVA, & FDCR are compared with derivatives calculated in a Fortran implementation using the complex-step method for accuracy. For N = 100, the ABACR method results in a term-by-term root-mean-square (rms) relative error of 10 −13 , while both FDSVA and FDCR result in an rms error of approximately 10 −12 . The rms relative error for all methods grows linearly with DoF on a log-log scale.
B. Runtime for Partials of ID/FD vs. RNEACR/FDCR as presented in [5] via Fortran implementation of both We consider an N link serial or branched kinematic tree with all revolute joints about their local z-axis.
RNEACR in body-coordinates [5, Algos. 2 & 3] (Table II) is used to calculate the partial derivatives of ID, and compared with the IDSVA (Table II) method. All the algorithms are written in Fortran 90 and implemented using the Intel Fortran compiler on a 3.07 GHz Intel Xeon processor. To calculate the average run time, each algorithm is run 10,000 times with randomized inputs for the state and control variables. Fig. 2 shows the comparison of the two methods. For N = 100, a speedup of 15× for IDSVA over RNEACR is found. Fig. 2 also shows the comparison of FDSVA with FDCR and ABACR (Table II)  Since SVA allows for coordinate-free expressions, recursive algorithms can be formulated in either bodycoordinates or ground-coordinates [19]. For the bodycoordinate algorithms, all the intermediate quantities are transformed between body coordinates using transformation matrices i X λ(i) , while for the ground-coordinates algorithms, the quantities can be left in the ground frame. A major advantage for the latter comes by avoiding the repeated transformation of the quantities between local body frames in the backward pass of the algorithm. This gain is achieved at the cost of expressing kinematic quantities (velocities 0 v i , accelerations 0 a i ), joint motion subspace matrices 0 S i , and inertia matrices 0 I i in the ground coordinate frame during the outward pass. Fig. 3 shows a comparison of IDSVA (see Table II) runtime in body coordinates and ground coordinates. Speedups for ground coordinate algorithms are between 1.3 to 1.9 for N = 2 to N = 500.

C. Comparing Runtime for C++ Partials of ID/FD with
Pinocchio Implementation Accompanying [5] IDSVA in ground coordinates is implemented in C++ within the Pinochhio framework. This strategy enables direct comparison with Pinocchio's original ID partial derivatives [6]. Fig. 4 shows the comparison of the two methods for serial and branched kinematic trees with some branching factor bf [19]. For a serial chain with N = 100, a speedup of 2× is found using the gcc-9.0  (Table II)

D. Comparing CPU Runtime for AZA vs. DMM in C++
with Pinocchio Implementation Accompanying [5] AZA is implemented within the Pinocchio [6] framework. Pinocchio's M −1 algorithm [25] is modified to include the AZA. In Fig. 7, the "crossover" N denotes the point below which the DMM performs better than the AZA. This point depends on the hardware and compiler optimization settings used to implement the algorithm, but for high N , the O(N 2 ) AZA is efficient because it avoids the expensive matrix-matrix product. Table III   gives the cross-over N values (for kinematic trees) above which AZA performs better than DMM for different compiler settings. The effect of this cross-over N can be seen in Fig. 6   calculate the inverse and forward dynamics partials for robots with multi-DoF joints and a floating base. Several algorithmic optimizations and Spatial Vector Algebra (SVA) identities are exploited to enable an efficient implementation. The method provides a 1.4× speedup for the TALOS humanoid model over the state-of-theart Pinocchio FD derivatives using the gcc compiler and a 1.2× speedup using the Clang-10 compiler. The reduction in runtime for partial derivatives enables faster optimization algorithms for both on-line and off-line applications. These improved timings can ultimately lead to better motion planning of legged and industrial robots.

A. Spatial Vector Algebra
A body k with spatial velocity k v k ∈ M 6 in the body frame is decomposed in its angular and linear components as: where k ω k ∈ R 3 is the angular velocity of the body in a coordinate frame fixed to the body, while k v k ∈ R 3 is the linear velocity of the body-fixed point at the origin of the body frame. Spatial vectors can also be expressed in the ground frame. For example, the spatial velocity of the body k in the ground frame is denoted as 0 v k . In this case, the linear velocity is associated with the body-fixed point on body k that is coincident with the origin of the ground frame. The net spatial force 0 f k ∈ F 6 defined in Eq. 36 on the body can be calculated from the spatial equation of motion (Eq.37): where 0 n k ∈ R 3 is the net moment on the body about the origin of the ground frame, 0 f k ∈ R 3 is the net linear force on body, 0 I k is the spatial inertia of the body k that maps motion vectors to force vectors, and 0 a k ∈ M 6 is the spatial acceleration of the body. The transformation matrix i X j is used to transform vectors in frame j to frame i is defined as: where i R j ∈ R 3×3 is the rotation matrix from frame j to frame i, p i/j ∈ R 3 is the Cartesian vector from origin of frame j to i, and 0 is the 3 × 3 zero matrix. p× is the 3D vector cross product on the elements of p, defined as: The spatial transformation matrix 0 X k can be used to obtain spatial velocity vector 0 v k from the vector k v k as: A spatial cross-product operator between two motion vectors (v,u), written as (v×)u, is given by (Eq. 41). This operation can be understood as providing the time rate of change of u, when u is moving with a spatial velocity v. A spatial cross-product between a motion and a force vector is written as (v× * )f , and defined in Eq. 42.
Hence, the three spatial vector cross-product operators defined above map between the motion and the force vector space as: [19] × :

B. Properties of Spatial Vectors
Assuming u, v, m, v 1 , v 2 ∈ M 6 , and f ∈ F 6 , many spatial vector properties [19] are utilized herein:

C. Multi-DoF joint Identities and Expressions
Identities are shown below for partial derivatives of common spatial quantities by perturbing a multi DoF joint position variable (q j ). These identities are used to then derive the partial derivatives of inverse dynamics for a rigid-body system with multi-DoF joints and a floating base.
otherwise Identities listed above are derived in detail as follows: J1. For directional derivative of S i along the p th freemode of the joint j, total derivative with respect to time is taken for the numerator and denominator as: Using the definition ofṠ i asṠ i = v i × S i , and the definition of v i : Only a single term remains for the case j i: J2. The directional derivative of the spatial velocity of a body i along the p th free-mode of the joint j, where j i is given as: Using J1, we get: By switching signs, Eq. 49 can also be written as: Eq 50 can be written for each p th mode of the joint j to get the partial derivative (of size 6 × n j , where n j is the number of DoF for joint j) of v i with respect to q j as: Using the definition of v l , Eq. 51 now is: Using Eq. 52, the partial derivative of v i × a with respect to q j for j i, where a ∈ M 6 is any fixed motion vector is given as: Similar to Eq. 53, the partial derivative of v i × * f , where f ∈ F 6 is any fixed force vector, is calculated as: J3. For any fixed force vector f , the partial derivative of S T i f with respect to q j for j i is calculated. Using J1: Eq. 56 can be written for each DoF of the joint j. Hence, the partial derivative with respect to q j is: J4. For any fixed motion vector a, the partial derivative of I i a is calculated with respect to q j for j i. The directional derivative of I i in the direction of p th free-mode of joint j [19] is: Multiplying a on both sides result in: Similar to Eq 56, Eq 59 can be written for each DoF of joint j collectively as: J5. The partial derivative of I i v i with respect to q j for j i is now calculated using the identities derived above. Using the product rule: For the first term in Eq. 61 on the RHS, v i is assumed to be a fixed motion vector. Hence, the identity J4 can be used. For the second term in Eq. 61, Eq. 52 is used as: Expanding terms: Upon cancellations and simplification, Eq. 63 becomes: whereΨ j is defined in Eq. 31. J6. Using identities derived above, and the definition of ξ i , the partial derivative of ξ i with respect to q j for j i is calculated as: where v J l is the joint velocity defined as: Using the product rule of differentiation: The partial derivative of the joint velocity v Ji with respect to q j for j i is calculated. Using the definition of v Ji (Eq. 68) we get: Using identity J1, we get ∂v Ji ∂q j,p = s j,p × S iqi (70) Collectively, for all DoF of joint j, Eq. 71 can be written as: Treating v J l as a fixed motion vector in the first term in Eq. 67 on the RHS, Eq. 52 is used. Eq. 53 is used for the second term.
Expanding terms, and using the definition ofΨ j (Eq. 31): Combining terms and simplifying, Using spatial vector cross property P1, definition of ξ i and simplifying, 76) Summing over the terms and simplifying, we get: J7. Using the definition of γ i , the directional derivative of γ i along the p th free-mode of the joint j for j i is calculated as: Using J1, we get: Eq. 79 can also be written as Collectively for all DoF of joint j , Eq. 80 is written as: Using the definition of γ l , Eq. 81 becomes: J8. The partial derivative of v i is calculated with respect toq j for j i. Using the definition of v i , we get: Summing over all indices l, all the terms vanish except the ones pertaining to the index j resulting in: J9. Using the definition of ξ i , the partial derivatives of ξ i with respect toq j for j i is: Using the product rule of differentiation, we get: Using J8, summing over the index l results in: Upon simplifying and using the definition ofΨ j (Eq. 31), we get: D. Partial Derivatives of ID w.r.t q: Derivations Partial Derivatives of [M (q)q] i : 1) Case when j i Using product rule of differentiation in Eq. 12, we get: Using identities J3, J4, and J7, we get: Upon cancellations and simplifications, the final expression is: 2) Case when j ≻ i For this case, identities J3, J4 and J7 are used as: Expanding: Upon cancellations and simplification, we get the following expression.
Partial Derivative of [C(q,q)q] i : 1) Case when j i Using product rule of differentiation in Eq. 13, we get: Using the property P3 to combine terms, Summing over the index k, we get: Partial Derivative of the Gravity Term:

1) Case when j i
For the case j i, the partial derivative of g i with respect to q j is: Using identities J3 and J4, we get: Cancelling terms, and summing over index k, we get: 2) Case when j ≻ i For the case j ≻ i, we follow the similar process and use identity J4 as: Expanding: Summing over the index k, we get: ∂g i ∂q j = S T i (I C j a 0 )× * S j + I C j (a 0 × S j ) (113) Re-arranging terms here UsingΨ k (Eq. 31), the final compact expression for ∂τ i ∂q j is: Re-arranging terms, we get: Using the definition of a i and f i (Eq. 10), we get: Using the definition ofΨ i (Eq. 31), we get the final expression for ∂τ j ∂q i as: F. Partial Derivatives of ID w.r.tq: Derivations 1) Case when j i Using product rule of differentiation in Eq. 13, we get: Using identities J8 and J9, we get: Simplifying and collecting terms, we get: Summing over the index k, we get: 2) Case when j ≻ i Similar to the previous case, product rule of differentiation in Eq. 13 is used: Using identities J8 and J9, we get: Simplifying and collecting terms, we get: Summing over the index k, we get: Flipping the indices i and j in Eq. 130 to get for the case j ≺ i: Hence, the partial derivatives of τ with respect toq are:


Recursive Bayesian Networks: Generalising and Unifying Probabilistic Context-Free Grammars and Dynamic Bayesian Networks

Probabilistic context-free grammars (PCFGs) and dynamic Bayesian networks (DBNs) are widely used sequence models with complementary strengths and limitations. While PCFGs allow for nested hierarchical dependencies (tree structures), their latent variables (non-terminal symbols) have to be discrete. In contrast, DBNs allow for continuous latent variables, but the dependencies are strictly sequential (chain structure). Therefore, neither can be applied if the latent variables are assumed to be continuous and also to have a nested hierarchical dependency structure. In this paper, we present Recursive Bayesian Networks (RBNs), which generalise and unify PCFGs and DBNs, combining their strengths and containing both as special cases. RBNs define a joint distribution over tree-structured Bayesian networks with discrete or continuous latent variables. The main challenge lies in performing joint inference over the exponential number of possible structures and the continuous variables. We provide two solutions: 1) For arbitrary RBNs, we generalise inside and outside probabilities from PCFGs to the mixed discrete-continuous case, which allows for maximum posterior estimates of the continuous latent variables via gradient descent, while marginalising over network structures. 2) For Gaussian RBNs, we additionally derive an analytic approximation, allowing for robust parameter optimisation and Bayesian inference. The capacity and diverse applications of RBNs are illustrated on two examples: In a quantitative evaluation on synthetic data, we demonstrate and discuss the advantage of RBNs for segmentation and tree induction from noisy sequences, compared to change point detection and hierarchical clustering. In an application to musical data, we approach the unsolved problem of hierarchical music analysis from the raw note level and compare our results to expert annotations.

Introduction
Long-term dependencies with a nested hierarchical structure are one of the major challenges in modelling sequential data. This type of dependencies is common in many domains, such as natural In contrast, DBNs are sequential models with a fixed set of random variables that reoccur at each time step [15,16]. The variables at each time step may be discrete or continuous, latent or observed, and may have an arbitrary non-cyclic dependency structure among each other, with additional links from the previous and to the next time slice. They comprise important model classes as special cases, such as hidden Markov models (HMMs) if there is only a single discrete latent variable or linear dynamical systems if all dependencies are linear Gaussians [17,18,16]. However, DBNs only allow for a fixed chain of Markov dependencies between time slices and cannot represent nested hierarchical structures.
In this paper, we present Recursive Bayesian Networks (RBNs), a novel class of probabilistic models that combines the strengths of PCFGs and DBNs by allowing for nested hierarchical dependencies in combination with arbitrary discrete or continuous random variables ( Figure 1). Our main contributions are as follows: 1. With RBNs, we provide a unified theoretical framework for a large class of important sequence models, including PCFGs and DBNs. 2. We generalise inside and outside probabilities from PCFGs to continuous latent variables, allowing for maximum posterior (MAP) inference in arbitrary RBNs. 3. For Gaussian RBNs, we derive an analytic approximation for the marginal likelihood and marginal posterior distribution, allowing for robust parameter optimisation and Bayesian inference. 4. We provide a quantitative evaluation on synthetic data and an application to the challenging task of hierarchical music analysis.

Related Work
PCFGs have a long tradition for modelling nested hierarchical dependencies in symbolic data with a variety of parsing algorithms for inferring the structure and variables' values [13,6,14]. Beyond their application to sequential data, PCFGs have been generalised to graph structures [19,20,21,22], which readily transfers to applications of RBNs. Latent vector grammars (LVeGs) [23] are an extension of latent variable grammars (LVGs) [24,25] with continuous latent states. As for RBNs, approximate parsing is possible in the Gaussian case. However, both LVGs and LVeGs are special cases of RBNs and do not draw the connection to graphical models. More recently, the availability of automatic differentiation libraries, such as PyTorch [26], has lead to a number of applications where gradients are propagated through the entire parsing process [27,28,29,30].
x z x x y z=N z=T Figure 2: RBN in Chomsky normal form. The non-terminal transition p N (x , x | x) and terminal transition p T (y | x) are grouped into an RBN cell with a structural distribution p S (z | x). We use gates [58] to describe structural distributions and extended factor graph notation [black squares; 59] for conditional joint distributions. Considering all possible ways how n observations can be generated by recursively applying the RBN cell produces an RBN chart, as shown in Figure 3.
Orange and red boxes, respectively, indicate the subsets X 0:3 and Y 0:3 of latent non-terminal and observed terminal variables generated from x 0:3 .
The process of parsing a PCFG or RBN can be formally rewritten as a sum-product network [SPN; 31,32,33]. Factor graph grammars [FGGs; 34] generalise PCFGs, case-factor diagrams [35] and SPNs by using a hyperedge replacement graph grammar [19] to describe a distribution over graph structures that is more general than that of RBNs (not only trees). However, none of the approaches addresses the problem of inference with continuous variables that we are facing in RBNs (exponentially many terms with exponentially many nested integrals).
A wide range of probabilistic and neural models operate with a fixed graphical structure and are loosely related to RBNs. Hidden tree Markov models [36,37] generalise HMMs from chain to fixed tree structures. They model data at each node as observations of a latent Markov process on the underlying tree, which is part of the input data. Additionally estimating the underlying tree structure has been addressed in [38,39,40]. Recursive neural tensor networks [41] use a PCFG for parsing a given sequence of symbols to obtain a tree structure, which is then fixed and used as the backbone for a neural network. More generally, there is a number of methods for inferring a fixed structure for graphical models [42,16,43,44], SPNs [45,46,47,48,32,49], or graph neural networks [50,51,52]. All these methods have in common that a fixed structure is either given or estimated but not treated in a probabilistic Bayesian manner.
Some approaches attempt a Bayesian treatment of the unknown structure of a graphical model or SPN via dynamic programming [53,54] or Markov chain Monte-Carlo sampling [55,56,57]. However, the structure is assumed to be independent of the latent variables (they only become dependent conditional on the data) and the latent variables cannot be used to control the structure, as it is the case in RBNs. The challenge of continuous variables also remains unsolved.

Recursive Bayesian Networks
RBNs are template-based graphical models that define a joint distribution over network structures and variables' values. The number of template variables is fixed, but the number of instantiated variables, their connectivity and values are governed by the joint distribution. As a rough analogy, RBNs can be thought of as DBNs that can not only be connected linearly to form a chain but also hierarchically to form a tree structure. Alternatively, they can be thought of as a PCFG in which each symbol is a (possibly continuous) random variable.

Definition
RBNs have three types of template variables: 1) latent non-terminal variables (discrete or continuous), 2) observed terminal variables (discrete or continuous), and 3) latent structural variables (always discrete). In the simplest case, illustrated in Figure 2, an RBN has one template variable of each type. Formally, an RBN is defined as follows: Definition 1 (Recursive Bayesian Network). An RBN is a tuple (X , Y, Z, T , S, p P ) with X : a set of latent non-terminal template variables (1) Y : a set of observed terminal template variables (2) Z : a set of latent structural template variables, paired up with the non-terminal variables (3) T : a set of transition distributions p(v 1 , . . . , v η | x) from a single non-terminal variable x ∈ X to a set of non-terminal and/or terminal variables v 1 , . . . , v η ∈ X ∪ Y (4) S : a set of structural distributions p(z | x), one for each non-terminal/structural pair (5) p P : a prior/start distribution for exactly one non-terminal variable.
The cardinality of a structural variable z ∈ Z corresponds to the number of possible transitions from the associated non-terminal variable x ∈ X ; η in (4) is called the arity of the transition.
Generating with an RBN is straightforward. We start by sampling the value of the first non-terminal variable x from the prior distribution p P (x) and then repeat the following steps until no unprocessed non-terminal variables are left: 1. sample the value of the associated structural variable from p(z | x) 2. choose a transition distribution p(v 1 , . . . , v η | x) based on the structural variable's value 3. sample the variables v 1 , . . . , v η from the transition distribution 4. for all newly generated non-terminal variables, go to step 1.
The major challenge and focus of this paper is to perform joint inference over the latent structure and non-terminal variables' values conditional on a given set of observations.
Chomsky Normal Form: In the simplest non-trivial case, an RBN has one latent non-terminal, one observed terminal, and one latent structural template variable, with one non-terminal transition of arity η = 2 and one terminal transition of arity η = 1, as illustrated in Figures 2 and 3. It is defined by four distributions In analogy to PCFGs, we call this the Chomsky normal form (CNF). Any RBN may be rewritten in CNF (see Appendix A.1 for details).
RBN Chart: During inference, we will make use of an RBN chart, similar to the parse chart for PCFGs [6]. Each non-terminal variable is associated to a layer in the chart. For discrete variables, they store the actual distributions, while for continuous variables they either hold the point estimate (for MAP inference) or the parameters of the approximate distributions (for inference in Gaussian RBN). Different instances of the same template variable are identified by a subscript indicating the span of data generated from them, which also corresponds to their position in the chart (see Figure 3). Sets of variables that are generated from a specific latent non-terminal variable x i:k are denoted by a bold capital letter with a corresponding subscript (X i:k , Y i:k , Z i:k ); omitting the subscript refers to all variables (X, Y, Z); for X and Z this also includes the root variables x 0:n and z 0:n , respectively. The subscripts are to be interpreted as time intervals, that is, Y i:i is empty, Y 0:1 = y 1 is the first observation, Y n−2:n = (y n−1 , y n ) are the last two observations etc.
Comparison to PCFGs: Any PCFG can be rewritten as an RBN in two different ways, which we call abstraction and expansion (see Appendix A.2 for details). Abstraction of a PCFG produces a discrete RBN with one latent non-terminal and one observed terminal variable. The resulting RBN is exactly equivalent to the original PCFG but describes the same relations in a more abstract and compact way. In contrast, expansion of a PCFG considers the symbols of the grammar as random variables in their own right, thereby endowing them with additional (possibly continuous) degrees of freedom. The resulting RBN is therefore more powerful than the original PCFG. A PCFG is abstracted to a discrete RBN by defining the start/prior, transition, and structural distributions (7-10) as where S is the grammar's start symbol, A, B, C are non-terminal symbols, b is a terminal symbol, X is any right-hand side of a rule, z=N and z=T indicate a non-terminal and terminal transition, respectively, W · → · is the weight of the corresponding rule, and rules that do not exist in the original PCFG are taken to have zero weight. In expansion, the PCFG is only used to define a "skeleton" for the RBN, while the specific random variables and the concrete transition distributions need to be additionally specified. This means that the resulting RBN model is more powerful than the original PCFG, as the symbols may, for instance, be expanded to continuous random variables.

Inference
The two main goals of inference in RBNs are to 1) train model parameters by maximising the marginal data likelihood and to 2) compute posterior distributions or maximum posterior (MAP) estimates of the network structure and non-terminal variables. In PCFGs, both is achieved by computing inside and outside probabilities [14], which will be the starting point for our generalisation to continuous variables.
Inside and Outside Probabilities: We define inside and outside probabilities, β and α, for RBNs in analogy to how they are defined for PCFGs, the only difference being that the variables may be continuous. We thus have where n is the length of the sequence and Y is fixed (and therefore omitted as argument on the left-hand side). That is, β(x i:k ) is the marginal likelihood of generating the sub-sequence Y i:k conditional on the respective non-terminal variable x i:k , while α(x i:k ) is the marginal likelihood of generating the two sub-sequences Y 0:i and Y k:n as well as the non-terminal variable x i:k . In both cases, β and α are functions of the corresponding non-terminal variable with the structure and the remaining variables being marginalised out. Based on the inside and outside probabilities, the marginal data likelihood and the marginal posterior distributions over non-terminal variables are p(Y) = β(x 0:n ) p P (x 0:n ) dx 0:n (17) and respectively. p(x i:k | Y) is an unnormalised probability distribution that specifies the probability of x i:k to exist via the normalisation constant p(x i:k | Y) dx i:k , while the normalised version corresponds to the marginal posterior distribution of x i:k for the case that it does exist.
Inside probabilities are recursively computed bottom-up. For an RBN in CNF we start with the base case (19) for single observations and then iterate (20) to the top of the RBN chart Outside probabilities are recursively computed top-down, while making use of the inside probabilities α(x 0:n ) = p P (x 0:n ) α(x j:k ) = p S (z j:l =N | x j:l ) p N (x j:k , x k:l | x j:l ) α(x j:l ) β(x k:l ) dx j:l dx k:l .
As for PCFGs, the two terms in (22) correspond to the possibility of x j:k being generated as the right or the left child, respectively. The main conceptual difference to PCFGs is that we treat the discrete structural part (marginalised out by the sums) separately from the potentially continuous variables (marginalised out by the integrals). For RBNs that are not in CNF, the equations have to be adapted accordingly (see Appendix A.3 for the general case).
Marginalisation: Computing the marginal data likelihood (17) and the marginal posterior distributions over non-terminal variables (18) requires to solve an exponential (w.r.t. the length n of the sequence) number of nested integrals in (19)(20)(21)(22), which is generally intractable. However, for the special case of Gaussian RBNs, we provide an adaptive closed-form approximation in Section 2.3. Moreover, marginalising only over the network structure for a fixed assignment of the non-terminal variables X is straight forward and allows for maximum posterior (MAP) inference in general RBNs.
Maximum Posterior Inference: For a fixed assignment of all non-terminal variables X, we can compute the joint marginal likelihood p(X, Y) over observed terminal and latent non-terminal variables by only marginalising over the structure. This follows the same principle as above but uses the modified joint inside and outside probabilities β i:k := p(X i:k , Y i:k | x i:k ) (23) and α j:k := p(X 0:j , Y 0:j , x j:k , X k:n , Y k:n ) , (24) where all variables are fixed (and therefore omitted as arguments on the left-hand side). Analogously, the joint marginal likelihood and the marginal posterior probability of x i:k to exist then are p(X, Y) = β 0:n p P (x 0:n ) (25) and where p i:k is the probability of x i:k to exist for this specific assignment of X. The corresponding equations for the recursion differ from (19)(20)(21)(22) only in that they do not integrate out the latent nonterminal variables (see Appendix A.3.3). As before, all computations can be efficiently performed via dynamic programming. Gradients w.r.t. the variables and/or parameters are readily obtained from libraries such as PyTorch [26]. Optimising the values of the latent non-terminal variables X via gradient descent yields maximum posterior (MAP) estimates, while the structure is marginalised out. MAP estimates for the structure (i.e. the best tree) conditional on an assignment for X can be computed (as for PCFGs) by replacing summation with maximisation [13,60].
There are two caveats: First, due to marginalising over multiple (exponentially many) network structures, p(X, Y) may be highly non-convex and optimising X via gradient descent is not guaranteed to find the global optimum. This is even the case for purely Gaussian RBNs, for which p(X, Y) is a mixture of Gaussians (one for each structure). Second, we can optimise X while marginalising out the structure and we can optimise the structure for a fixed assignment of X. However, successively optimising X and the structure is not equivalent to jointly optimising both and the maximum of p(X, Y) may be unrelated to the maximum of the best structure (also see Figure 4). This means that generally, exact joint MAP inference over the latent variables and the structure is hard. For Gaussian RBNs, we provide an approximate solution below.

Gaussian RBNs
In a Gaussian RBN (GRBN), the prior, non-terminal, and terminal distributions are linear Gaussians and the termination probability (structural distribution) is constant [termination/structural] For clarity, we will show all derivations for GRBNs in this basic form. For our evaluations and the application to music, we use a slightly extended version that includes linear transformations, mixtures of Gaussians, and multi-terminal transitions (Section 2.3.1). The derivations do not fundamentally change for the extended case (see Appendix A.4). In Appendix B, we show all calculations on a simple example.
Adaptive Approximation: If the structure of a GRBN was fixed, all variables would be jointly Gaussian distributed as in a conventional Gaussian Bayesian network [16]. However, due to the mixture max approximation components Figure 4: Three Gaussians components, the resulting mixture (green), maximum (orange), and moment-matching single Gaussian approximation (blue). Note that the maximum of the mixture (), the best component (), and the approximation () may be unrelated.
x λ W x x τ y n n z=N z=T z Figure 5: Graphical model of the Gaussian RBN for modelling music. The additional transposition variable τ is marginalised out during inference; the number of jointly generated observations n is uniquely determined by the location in the parse chart.
unknown structure, we effectively have a mixture of exponentially many Gaussians, one for each possible structure. While in principle all integrals can be solved analytically, the exponential growth makes exact inference intractable. Therefore, our goal is to derive a parsing strategy that retains tractability by adaptively applying local approximations to the Gaussian mixtures occurring in each recursion step. We will here focus on the simplest case of approximating the mixtures with a single Gaussian (illustrated in Figure 4, details in Appendix A.4.2), which can be efficiently computed in closed form [18,61]. The inside and outside probabilities are thus represented by a simple Gaussian and this form is reestablished in each iteration by approximating the occurring mixtures. Consequently, the marginal posterior distributions over latent variables (18) are also simple Gaussians and the marginal data likelihood (17) can be computed in closed form. This approximation scheme can be extended and refined by using existing methods for approximating each Gaussian mixture by one with fewer components [62,63].
Marginalisation: In (20) and (22), we have to integrate over products of Gaussian distributions to marginalise out the latent variables. To solve these integrals, we make use of the fact that the product of two Gaussians over a variable x can be rewritten as [see e.g. 64] withΣ Hence, when integrating over x, only the first term on the rhs. of (33) remains. A detailed step-by-step derivation of all results can be found in Appendix A.4.1. With the latent variables being marginalised out, (20) and (22) become simple mixtures of Gaussians that can be easily approximated to retain the simple analytic form of the inside and outside probabilities.
Tree Induction: As described above, exact joint MAP inference over the continuous latent variables and the structure is generally intractable. Moreover, the maximum of the approximate posterior does not necessarily coincide with the maximum of the exact posterior or that of a particular structure (see Figure 4). Thus, first optimising X (based on the approximation) and then estimating the structure (conditional on the picked value of X) may lead to arbitrarily bad results for tree induction. Therefore, we leverage the adaptive character of our approximation scheme to compute local structure estimates in each step, before loosing relevant information due to further approximations. Specifically, during the bottom-up pass for computing inside probabilities, all structures are scored by the maximum of their marginal likelihood, based on its current approximation (31). The best overall structure is then selected (as usual) in a top-down pass (see Appendix A.4.3 and our example in Appendix B).

Gaussian RBNs for Music
For the application to music, we slightly extend the basic GRBN discussed so far by introducing transpositions and multi-terminal transitions (changes in the equations highlighted in blue). The corresponding graphical model of the RBN cell is shown in Figure 5. Furthermore, we describe how GRBNs can be applied to categorical data.
Transpositions: A transposition rotates the dimensions of the latent variable by a number of steps τ before generating the child. This is achieved by multiplying with an orthonormal transposition matrix T τ that corresponds to the identity matrix with cyclicly rearranged columns. For the prior distribution, we assume a uniform weighting of all possible transpositions where D is the dimensionality of the data (D = 12 for music in 12-tone equal temperament). For the non-terminal transitions, the probability for a specific transposition is determined by the weight parameter W Note that transpositions are only applied to the left child, because Western classical music is thought to be fundamentally goal directed [2,8,65]. This means that the character of a section is largely determined by how it ends (the right child), which should also be reflected in the value of the parent node. In contrast, the role of the left child is to harmonically prepare the ending (or prepare a preparation to the ending etc). We therefore allow for arbitrary transpositions in the left child and we will see below that our model indeed captures the most important type of preparation in Western classical music: the cadential dominant-tonic progression.

Multi-Terminal Transitions:
A multi-terminal transition generates multiple observed variables from a single latent variable. The variables are generated i.i.d. and their number is governed by a Poisson distribution with rate parameter λ Multi-terminal transitions do not conform to the CNF assumed so far and we need to add the term [for GRBNs, see (30)] (39) to (20) in order to account for the possibility to terminate from a higher-level variable. For k = i + 1, this term becomes the base case (19) of an RBN in CNF.
Multi-terminal transitions account for the situation where changes in the hierarchical structure occur at a lower rate than the time series is sampled. In between the structural changes, the data is assumed to be generated from the same model, which could also be more elaborate than i.i.d. samples, as long as the relevant model parameters are captured by the RBN's latent variables.
Categorical Data: The observed variables of a GRBN are unconstrained real-valued, which poses a problem if the data are categorical. This situation is comparable to using Gaussian processes (GPs) [66] for classification and can be approached with similar methods. In our application to musical data, we observe one or more notes being played at any particular time and normalise these counts to obtain observations that correspond to the parameter of a categorical distribution. The natural likelihood function for this type of observations is a Dirichlet distribution. Therefore, we adapt the approach suggested in [67] for GPs, who assume a Dirichlet likelihood, which is then approximated by a Gaussian likelihood in log-space. Since an observation from a Dirichlet distribution corresponds to a normalised sample from independent Gamma distributions, each Gamma distribution can be separately approximated by a log-normal distribution, which results in a diagonal covariance matrix for the Gaussian likelihood in log-space. Matching the first and second moment yields [67] y (l) where 0 < y (l) j < 1 is the l th element (normalised count) of the j th observation, y j is the corresponding mean of the approximate Gaussian likelihood in log-space, and Σ (j) ll is the l th element on the diagonal of the covariance matrix for the j th observation. We thus have to replace y j and Σ (j) for y j and Σ T in (37).

Experiments
We performed a quantitative evaluation on synthetic data and applied our model to hierarchical music analysis of Bach preludes. We show that RBNs are superior to change point detection (CPD) and hierarchical clustering (HC) for tree induction and our method is able to infer fundamental harmonic principles of Western classical music. Experiments were run on a 3.6 GHz Quad-Core Intel Core i7 processor with 32GB RAM. The model parameters were trained via gradient descent on the (approximate) marginal neg-log likelihood.

Quantitative Evaluation on Tree Induction
We performed a quantitative evaluation on synthetic data for the task of segmenting a noisy time series and inferring the underlying tree. For comparison, we used the best-performing change point detection (CPD) method from the ruptures library [72] for segmenting the time series, combined with bottom-up hierarchical clustering (HC) for inferring the tree structure ("HC/CPD"). For details of the methodology, see Appendix C.1.
The evaluation results in Figure 6(a) show that the RBN tree estimates (Section 2.3) consistently outperform the one from HC/CPD, in terms of both precision and recall (and thus also in F1 measure). The marginal node probabilities show an interesting performance pattern. They excel in terms of precision, which means that a node with high marginal probability is very likely to actually exist in the tree (low false-positive rate). However, they severely underestimate the overall node probabilities, which leads to recall falling far below the baseline. This means that a node with low marginal probability may in fact occur in the tree (high false-negative rate).
We think that the poor recall measure of the marginal probabilities is primarily due to (and the downside of) a fully Bayesian treatment that quantifies uncertainty. Even if the marginal probabilities have a maximum at the correct node location, probability mass will still spread around it and be allocated to a number of less probable locations. While this is the desired behaviour of a Bayesian method, it inevitably results in a lower recall value. The high precision value confirms that uncertainty is adequately quantified and not underestimated. That being said, the marginal probabilities provide an exceptionally rich basis for qualitative analyses. For instance, all ground-truth nodes are located at local maxima of the marginal probabilities and we can read off a number of other potential node locations, which essentially trace out the grid defined by the piece-wise constant segments (see Figure 8 in Appendix C.1).

Hierarchical Music Analysis
Harmonies in Western classical music exhibit a nested hierarchical structure that can be modeled by PCFGs operating on abstract chord symbols [8,73,10,3]. While these grammars can be applied to expert annotations of a musical score, hierarchical music analysis from the raw note level is an unsolved problem. We trained a GRBN (Section 2.3.1) on the 24 major preludes of Johann Sebastian Bach's "Wohltemperiertes Klavier I & II" (see Appendix C.2 for technical details and complete results).
Our first major finding is that the prior mean, shown in Figure 6(b), corresponds to a major pitch profile (as could be expected from the training data) and is in excellent agreement with recent Bayesian estimates from the literature [68]. The fact that the major profile appears in the prior (i.e. as the continuous equivalent of a grammar's start symbol) shows that our model picks up fundamentally important structures from the musical data. Our second finding is that only two transpositions have non-zero weights: the identity with a weight of 78% and the fifth scale degree (7 semitones) with a weight of 22%. This corresponds to the left child being generated as the dominant of the parent and realises the most important harmonic preparation in Western classical music: the cadential dominant-tonic relation. A closer inspection of the expert analysis ( Figure 10 in Appendix C.2) reveals that when considering the possible surface patterns (raw notes) of the labeled chords, most non-identity transitions can indeed be explained as (noisy) fifth transpositions. The strong weight of fifth transpositions in our model is a highly non-trivial empirical confirmation of the established music theoretical insight that Baroque music is fundamentally driven by dominant-tonic relations. While the estimated tree in Figure 6(c) fails to reproduce the large-scale structure of the expert analysis (e.g. the separation into two main parts), it accurately captures the measure-wise harmonic changes on the bottom level.
On the one hand, we see considerable room for improvement by integrating more advanced concepts, such as different modes (major/minor), diatonic in addition to chromatic transposition, or balancing of trees. On the other hand, our model was able to capture fundamental properties of Western classical music based on only 24 pieces. We therefore think that Gaussian RBNs are a highly promising approach for hierarchical music analysis from the raw note level, which should be further investigated.

Conclusion
We introduced Recursive Bayesian Networks (RBNs), a novel class of probabilistic models that unifies the strengths of probabilistic context-free grammars (PCFGs) and dynamic Bayesian networks (DBNs), generalising both model classes. We defined RBNs as a joint distribution over tree-structured Bayesian networks and their (discrete or continuous) variables and described how to perform inference over both the model structure and the variables by leveraging parsing methods for PCFGs. The provided formalisation connects with the methods for formal grammar as well as with the versatile notation for graphical models. On two data sets, we demonstrated the potential of RBNs for modelling nested hierarchical dependencies in real-valued time series and musical data. The class of RBNs represents a substantial contribution to the machine learning toolkit by unifying two of the most important approaches for modelling sequential data and bears a large potential for further development and applications.

Appendix
Recursive Bayesian Networks: Generalising and Unifying Probabilistic Context-Free Grammars and Dynamic Bayesian Networks Transforming an RBN into CNF is done analogously to the procedure for PCFGs. We assume the original RBN does not contain any epsilon productions, that is, a non-terminal variable always produces one or more other non-terminal and/or terminal variables.

1) Eliminate terminal variables from mixed transitions: This is done by introducing intermediate non-terminal variables. For each transition
with non-terminal variables x (1) , x (2) , . . . ∈ X and terminal variables y (1) , y (2) , . . . ∈ Y, we introduce new non-terminal variables x y (1) , x y (2) , . . . and replace the transition by where the new non-terminals x y (i) replace the original terminals y (i) . We then add new deterministic transitions . . . that convert each non-terminal to its equivalent terminal variable. For the newly added non-terminals, there is only a single transition and hence a degenerate structural variable that can only take a single value.
2) Eliminate more than two latent non-terminal variables: This is done by introducing new nonterminals that capture combinations of multiple old non-terminals. Below, we show how the number of non-terminals can be reduced by one. Applying this procedure repeatedly allows for reducing the number of non-terminals from an arbitrary number down to two, as required for CNF. A transition that generates n non-terminals x (1) , . . . , x (n) is rewritten as where we introduced the new non-terminal variable x = (x (1) , x (2) , . . . , x (n−1) ) that stores all the information from the first n − 1 original non-terminals. The actual "work" is done by p(x , x (n) | x), which is the equivalent of the original n-fold transition. p(x (1) , x (2) , . . . , x (n−1) | x ) is a deterministic transition that just "unpacks" the information stored in x . Repeating this procedure to come to only pairwise transitions corresponds to a chain of these deterministic "unpacking" operations. As above, the newly added non-terminals have only a single possible transition.
3) Eliminate unary cycles: Unary cycles p cycle (x | x), where x and x are the same non-terminal template variable (x ≡ x ), are first transformed into unary transitions to a new non-terminal variable and then eliminated as described below. We define a new non-terminal variablex = (x , n), where n > 0 represents the number of steps taken in the cycle before exiting it and x is the value at the moment of exiting it. The transition distribution tox is where in the recursive case, the variables of all intermediate steps are successively marginalised out. In practical applications, if p(z=cycle | x ) < 1, the probability of remaining in the cycle decays exponentially and the recursion can be truncated after a number of steps. If on the other hand p(z=cycle | x ) ≈ 1 so that truncating is not possible, one can work with the stationary distribution of the resulting Markov chain (i.e. the Markov chain with transition distribution p cycle (x | x)).
The structural probability to take a transition from x tox is p(z=cycle | x), i.e. the probability of entering the cycle in the first place. The RBN cell ofx is identical to that of x, except for the transition into the cycle, which is eliminated (the structural distribution thus has to be renormalised for the remaining transitions). The transitions use only the x -component ofx, ignoring the n-component. In this way, we have expresses the state after an arbitrary number of steps in the unary cycle as a distinct value of the new non-terminal variablex.
The intermediate variable x and its RBN cell is eliminated if it was only reachable via x. The new transitions p 1 , . . . , p n are merged into the cell of x, while the original transition p unary to x is removed. This requires redefining the structural distribution p(z | x) such that the probability mass p(z=unary | x) that was formerly assigned to p unary is now split among the new transitions p 1 to p n according to the structural distribution p(z | x ) of x . Specifically, for a new transition p i , we define A.2 Relation to PCFGs As described in Section 2.1, a PCFG can be rewritten as an RBN by abstraction or expansion, where abstraction produces an equivalent RBN that describes the same relations in a more abstract and compact way, while expansion produces a more general RBN using the original PCFG as a skeleton.
We describe the two procedures in detail below and use the following definition of a PCFG: where A, B, C ∈ N are non-terminal symbols of the PCFG, b ∈ T is a terminal symbol, X ∈ N 2 ∪T is any right-hand side of a rule, z=N and z=T indicate a non-terminal and terminal transition, respectively, W A→X is the weight of the corresponding PCFG rule, and rules that do not exist in the original PCFG are taken to have zero weight.
To show equivalence, we need to prove that the transition probabilities from a given non-terminal symbol A are the same in the original PCFG and the new RBN.
Proof. In the RBN, the probability for a non-terminal transition A → B C is and that for a terminal transition A → b is which matches the corresponding probabilities in the PCFG, gained by normalising the respective weights.
Conversely, any discrete RBN can be rewritten as a PCFG. Theorem 2. A discrete RBN with n latent non-terminal template variables x 1 , . . . , x n , m observed terminal template variables y 1 , . . . , y m , and a prior p P (x 1 ) over x 1 can be rewritten as a PCFG with where · ⊕ · concatenates the value ranges of the respective variables, X ∈ x i denote that the value X is in the value range of the RBN variable x i , and the second case in (60) requires there be a transition p i (x 1 , . . . , x k | x i ) such that A ∈ x i and X ∈ x 1 ⊕ · · · ⊕ x k .

A.2.2 Expansion of a PCFG
Expansion of a PCFG to an RBN uses the PCFG as a "skeleton" to define the number of template variables and the structural transitions. The domains and transitions for the variables need to be added, which results in an RBN that is more powerful than the original PCFG. Specifically, we have for the sets of latent non-terminal and observed terminal template variables and for the structural transitions. Additional, we have to define the domain for each of the non-terminal and terminal variables in X and Y, and for each rule A → X 1 X 2 . . . from the original PCFG, we have to define a concrete transition distribution . . ∈ X ∪ Y are non-terminal or terminal variables in the RBN, respectively, depending on whether X 1 , X 2 , . . . ∈ N ∪ T are non-terminal or terminal symbols in the PCFG).
Expansion of a PCFG into an RBN seems appealing if a simple PCFG can be used to describe the type of variables (as opposed to their values) as well as the structure of the generative process. The actual transitions on the variables' values may then take place on a sub-symbolic/continuous level, which cannot be described by a PCFG.

A.3.1 Inside Probabilities
The inside probability is the probability of generating the observed terminal variables Y i:k from the latent non-terminal variable x i:k . This means that we need to marginalise over all possible paths of generation. Transitions may directly generate observed variables, but they may also generate lower-level non-terminals, in which case we have to recurse using the respective inside probabilities from those variables.
Let T x ⊆ T be the set of possible transitions from the latent non-terminal template variable x ∈ X (of which x i:k is one specific instantiation), with p(z i:k =τ | x i:k ) being the probability for the transition τ ∈ T x to be selected. This constitutes the first sum in (64) below, which marginalises over the transitions.
The transition τ generates η new non-terminal and/or terminal variables, where η is the arity of τ . These may be located at different positions in the parse chart, depending on which part of the observed variables Y i:k is generated from them. That is, the variables' locations in the parse chart are not known during generation and are determined in hindsight once all observed variables are generated; thus, they are known for parsing. We denote the respective splitting points by j 1 , . . . , j η−1 (they have to fulfill the condition i < j 1 < . . . < j η−1 < k) and the respective variables by v i:j1 , . . . , v jη−1:k ∈ X ∪ Y. The second multi-sum in (64) aggregates the probabilities of the different splitting possibilities, that is, of all valid assignments of j 1 , . . . , j η−1 (η − 1 degrees of freedom). For instance, a transition of arity η = 2 has one free splitting point j 1 to sum over.
Some of the generated variables may be observed/terminal variables, for which nothing more needs to be done as they directly constitute the respective part of Y i:k . For the subset of non-terminal variables, which we denote by {v j:j ∈ X }, we need to insert their respective inside probabilities and marginalise them out. This constitutes the product and multi-integral in (64).
The general form of the inside probabilities then is The concrete RBNs considered in the paper have only two transitions, one non-terminal transition of arity two and one terminal transition of arity one (for CNF) or more (for the extended GRBNs used in the quantitative evaluation and for modelling music). For non-terminal transition of arity two, the multi-sum in (64) reduces to a single sum and the multi-integral to a double integral, which gives us (20). For the terminal transition, (64) simplifies to (19) or the extended version (38), respectively.

A.3.2 Outside Probabilities
The outside probability is the joint probability of generating the latent non-terminal variable x j:j as well as the prefix and suffix of observed terminal variables, Y 0:j and Y j :n , respectively. For this, we now have to consider all possible ways how x j:j as well as the prefix and suffix could have been generated from a parent non-terminalx (x andx may correspond to the same template variable or to two different ones).
Let T −1 x ⊆ T denote the set of transitions that include x as one of the generated variables. Importantly, if x appears multiple times in the generated variables of a particular transition, these different options of generating x are represented as multiple distinct entries in T −1 x , one for each occurrence. The first sum in (66) runs over these different possibilities of generating x.
For a transition τ ∈ T −1 x of arity η, let j 0 , . . . , j η be the splitting points, including the start and end point j 0 and j η of the parent variablex j0:jη , which have to fulfill the condition 0 ≤ j 0 < . . . < j η ≤ n (where n is the length of the sequence). One pair of adjacent splitting points (j m , j m+1 ) corresponds to the occurrence of x j:j , where m is the position (starting at zero) at which x appears in the generated variables of the particular transition τ . We therefore have the additional constraints j m = j and j m+1 = j , resulting in η − 1 remaining free indices to sum over (as for the inside probabilities above). This corresponds to the second multi-sum in (66).
The set of non-terminal variables generated from the parentx j0:jη , excluding x j:j , is denoted by {v l:l ∈ X } \ x j:j . Together with the directly generated terminal variables, these generate part of the prefix and suffix, Y j0:j and Y j :jη . The remaining prefix and suffix, Y 0:j0 and Y jη:n , are generated from the parent variablex j0:jη . For the parent, we recurse via its outside probability α(x j0:jη ), while for the newly generated non-terminal variables (except x j:j ), we have to use the respective inside probability β(v l:l ) in (66). Additionally, we have to marginalise out the parent (first integral) and the newly generated non-terminal variables (second multi-integral). For a non-terminal transition of arity two, as we have it in the paper, the multi-sum in (66) reduces to a single sum and {v l:l ∈ X } \ x j:j contains only a single non-terminal, the second child. Importantly, T −1 x has two elements, one for x j:j being generated as the right child and one for it being generated as the left child, which gives us (22).

A.3.3 Joint Inside and Outside Probabilities
The joint inside and outside probabilities (23) and (24) for an RBN in CNF are computed analogously to (19)(20)(21)(22) for the normal inside and outside probabilities, that is, This differs from (19)(20)(21)(22) only by dropping the integrals and dependencies on the non-terminal variables (as their values are now fixed). Joint inside and outside probabilities for the general case are obtained from (64) and (66) analogously, i.e. again by dropping the integrals and dependencies on the non-terminal variables.

A.4 Gaussian RBNs
In the following, we present derivations for the extended case of GRBNs, described in Section 2.3.1, which includes linear transformations T for the left child. For this, we will make use of the fact that a normal distribution over a transformed variable T x can be rewritten as where ||T || is the absolute value of the determinant of T and in (72) we made use of the fact that in our case, the transformation matrices are orthonormal, so that T −1 = T and ||T || = 1.
Note that for an implementation, some of the results should be rewritten in order to minimise the number of matrix inverses that need to be taken. In particular, the identity is useful for the implementation, but we omit it in our derivation for clarity.

A.4.1 Marginalisation
For the inside probability β(x i:k ), the integral in (20) is N (x i:j ; µ where in (74) we inserted (31) and (36); in (75) we used (33) twice to rewrite the pairwise products of Gaussians over x i:j and x j:k and marginalise them out; in (76) we used (72) to rewrite the transformation; and in (77) we used (33) a third time to rewrite the resulting product as a single Gaussian over x i:k .
For the outside probability α(x j:k ), the integrals in (22) for x j:k being generated as the right child are with and where in (81) we took the constant termination probability (30) out of the integral and dropped it; in (82) we inserted (31), (32) and (36); in (83) we applied (33) to marginalise out x i:j ; in (84) we used (72) to rewrite the transformation; in (85) and (86) we used (33) twice to marginalise out x i:k ; and in (87) we rewrote the final result using (91-93). Due to the asymmetric terms in the outside probabilities, the result is somewhat more complex than for the inside probabilities.
Analogously, the integrals in (22) for x j:k being generated as the left child are p N (x j:k , x k:l | x j:l ) α(x j:l ) β(x k:l ) dx j:l dx k:l (94)

A.4.2 Approximation
A Gaussian mixture distribution p(x) with normalised mixture weights c i , means µ i , and covariance matrices Σ i can be approximated with a single Gaussian as The approximation p(x) matches the first and second moments of p(x) and minimises the Kullback- [61,18]. This direction of the KLD is the one used e.g. in expectation propagation, not the one used in e.g. variational methods [18]. That means, p(x) will adequately represent the support and uncertainty of p(x) (e.g. it will be non-zero wherever p(x) is non-zero). On the other hand, a value of x may have a high probability in p(x) even though in p(x) it has not (also see Figure 4).

A.4.3 Tree Induction
Exact joint optimisation of the structure and the continuous latent variables is intractable. We therefore choose the best tree for a GRBN based the maximum of the (approximated) inside probability. Inserting (30) and (77) x 0:4 x 0:3 x 3:4 x 0:1 x 1:3 x 1:2 x 2:3 y 1 y 2 y 3 y 4 Figure 7: Parse chart for a sequence of length n = 4 with best tree estimate (see text for details).
which is maximised by taking the mode of the Gaussian and maximising over j and (if using If we have multi-terminal transitions (or more generally other possible transitions), we also have to maximise over the different possible transitions. For each non-terminal variable, we compute and store the best choice during bottom-up computations of the inside probabilities. Afterwards, we can construct the best tree by starting at the root node and recursively picking the best structure top-down.

B Example
In this section, we present the complete calculations for the inside probabilities, the tree estimate, and the marginal likelihood for a basic GRBN (no transpositions or multi-terminal transitions) on a simple example sequence of length n = 4 (also see Figure 7). We assume parameters in (27)(28)(29)(30) and a scalar sequence Y = (y 1 , y 2 , y 3 , y 4 ) = (0, 1, 2, 0) .

B.1 Preliminaries
Inside probabilities are approximated with a single Gaussian specified by c i:k , which are the relevant quantities to be computed. At the bottom level, we use (19) for the base case and insert (29) and (30) to obtain where we can directly read off c where we already inserted Σ NL = Σ NR = 1.
If the sum in (113) has only a single term, we immediately get If there is more than one term in the sum in (113), this means that there are multiple splitting options that are marginalised out and we therefore need to do two things.
First, we have to identify the best splitting option to be able to compute the tree estimate. This is done by using (109) and comparing the values of where Σ Second, we have to approximate the resulting mixture with a single Gaussian using (107), where the mixture weights have to be normalised. For the univariate/scalar case considered here, we then get Finally, the marginal likelihood (17) where we have used (33) to rewrite the product of Gaussians and inserted µ P = 0 and Σ P = 1.
The normal distribution is defined as where the second line is for the scalar case.

B.2 Calculations
We start with the inside probabilities at the bottom level for the latent variables x 0:1 , x 1:2 , x 2:3 , x 3:4 and from (112) we read off (without any approximations) Next, we compute the inside probabilities on the first level for the variables x 0:2 , x 1:3 , x 2:4 . The only possible splitting option is for j = i + 1 and from (117) we get (again without approximation) Turning to the values for x 0:3 and x 1:4 , we now have two terms in the sum in (113), which means that we need to evaluate the best split and approximate the mixture. The corresponding parameters of the mixture are given by (114-116) as that for x 0:3 both splits are equally well and for x 0:3 the split x 1:4 → (x 1:3 , x 3:4 ) at j = 3 is better. This is intuitively clear, since generating (y 2 , y 3 ) = (1, 2) from the same non-terminal variable x 1:3 = 1.5 is more likely than generating (y 3 , y 4 ) = (2, 0) from x 2:4 = 1, given that in both cases the values are generated from a Gaussian with variance 1.
We approximate the mixtures with a single Gaussian with parameters given by (119-121) as and we see that the split x 0:4 → (x 0:3 , x 3:4 ) for j = 3 is the best one. Intuitively, this makes sense because it splits between y 3 and y 4 , which is the biggest step. Not we have only minor differences between the split options, because for simplicity we have chosen our variance parameters with a value of 1, which is relatively large compared to the spread of the values. Choosing smaller variances would result in more prominent splitting preferences.
We can now construct the full tree by also picking the best split for x 0:3 , which is a tie between splitting at j = 1 and j = 2, so we can choose either one (in practice one might consider random tie breaking to avoid biases due to variable order). The resulting tree is shown in Figure 7.
The parameters for the inside probability of x 0:4 , given by approximating the three Gaussian mixture components, are

C.1 Details for Quantitative Evaluation
We performed a quantitative evaluation on synthetic data for the task of segmenting a noisy time series and inferring the underlying tree. We used the Gaussian RBN for music (Section 2.3.1) with  some simplifications: 1) data were continuous and not categorical, 2) they had only three dimensions instead of twelve, 3) the prior distribution did not have any transpositions, 4) the left child could only be transposed by zero or one step. The data were sampled from this model, using zero prior mean, Σ p = 1, Σ NL = Σ NR = 1 · 0.1 2 , λ = 5, equal weights W 0 = W 1 = 0.5 for transposition by zero or one, and Σ T = 1 · noise 2 with different noise levels {0.01, 0.05, 0.1, 0.15, 0.2, 0.25}. We used a terminal probability of p term = 0.6 for sampling and rejected any sequences with a length outside the range of 50-55. An example of the data is shown in Figure 8.
For comparison, we used the best-performing change point detection (CPD) method from the ruptures library [72] for segmenting the time series, combined with bottom-up hierarchical clustering (HC) for inferring the tree structure ("HC/CPD"). For each noise level, we selected the CPD method and parameters with best F1 score based on the ground-truth segments of 100 training sequences (see Figure 9). In HC, pairs of adjacent segments with the smallest Euclidean/L2 distance between their mean values were successively combined to construct the tree. For the RBN, all parameters were trained from scratch by minimising the marginal likelihood of the observations of only 10 training sequences (i.e. no ground-truth information and less training data than for the baseline was used), separately for each noise level.
The models were evaluated on 500 test sequences by computing their precision and recall w.r.t. the ground-truth trees. Each possible node was treated as a separate binary classification task and the results reflect the number of correctly predicted nodes. Note that due to the strongly unbalanced class distribution (many more possible node locations than actual nodes in the tree) precision and recall or the combined F1 score are the appropriate performance metrics (as opposed to e.g. accuracy). For the RBN they were computed in two different ways: 1) based on the best-tree estimate ("RBN max") and 2) based on the marginal node probabilities (18) For the single-tree estimates (baseline model and best-tree estimate from RBNs) we compared the ground-truth and estimated tree node-by-node to count correctly predicted nodes (TP), nodes that are in the prediction but not the ground-truth (FP), and nodes that are in the ground-truth but not the prediction (FN). For the marginal node probabilities, we computed the corresponding rates by counting all nodes in the ground-truth tree (TP+FN), summing the marginal probabilities over all

C.2 Hierarchical Music Analysis
The scores were pre-processed by computing pitch-class distributions (PCDs), as used for the identification of musical keys [74,75,76], using the pitchscapes library [71]. We used a resolution of 70 equally spaced time slices per piece, resulting in sequences of 12-dimensional categorical distributions. The tree shown in Figure 6(c) for Johann Sebastian Bach's Prelude No. 1 in C major, BWV 846, corresponds to a harmonic expert analysis performed by the authors. Figure 10 shows the annotated tree with additional chord labels, which are provided in a simplified notation commonly used in Jazz lead sheets to be more accessible to a broad audience. In Table 1, we list the results for all 24 preludes. For a better interpretation of the model and the presented results, there are two relevant points to consider.

C.2.1 Chromatic versus Diatonic Transposition
It is interesting to look in more detail at what musical aspects the model can or cannot represent. In a nutshell, it can represent chromatic transposition but cannot represent diatonic transposition, which has a number of consequences, as described in the following.
The transpositions of the left child perform a cyclic rotation of the corresponding probabilities in the pitch-class distribution represented by the latent variable, which corresponds to a chromatic transposition. This determines not only which pitch classes have a significant probability to occur (the in-scale tones) but also the specific weights. For instance, the tonic and fifth scale degree typically have the highest weights. A transposition by 5 or 7 semitones from a current major key (say C major) thus corresponds to a modulation to the sub-dominant (F major) or dominant (G major) key, respectively. This includes adaptation of the fourth and seventh scale degree of the target key, respectively (B→B for F major; F→F for G major), as well as the correct assignment of strong weights to the tonic and fifth scale degree.
However, diatonic transposition cannot be represented in this way. For instance, to represent a modulation from C major to A minor, the model has two options that are both far from optimal. 1) It can choose not to apply a chromatic transposition, which ensures that all in-scale tones are correctly represented (i.e. they have significant weight). This, however, means that the relative weights are not appropriate for A minor. In particular, the strong weights on the tonic and fifth scale degree are not present and, instead, the third and seventh scale degree (C and G, the former tonic and fifth scale degree) have disproportionally strong weight. Correcting these weights has to occur through the Gaussian transitions, which can only be explained with a relatively high transition variance.
2) The second option would be to perform a chromatic transposition by 9 semitones, which ensures that the strongest weights remain on the tonic and fifth scale degree of the new key. However, three out-of-scale tones (C , F , G ) now have a high weight, while the respective in-scale tones do not. Again, this has to be corrected for by the Gaussian transition noise at a potentially even higher cost than in the first case. This is a highly plausible explanation for why we only see non-zero weights for the identity and (chromatic) transposition by a fifth in our experiments. Any diatonic modulations are best explained by reweighting using via Gaussian transition noise without a transposition, rather than by a chromatic transposition, which would require an even stronger reweighting (except for modulation to the sub-dominant and dominant key, which can be appropriately explained by a chromatic transposition).

C.2.2 Chord Labels
It is important to note that the chord labels in the expert annotation convey significantly more information than just what pitch classes can be expected to occur in the respective section. For example, the very same pitch-class distribution of G-C-E could amongst others be labeled as a C major chord in second inversion, a G major chord with 64-suspension (Gsus64), or an A minor seventh chord with omitted root, which might be easily confused by a musically untrained annotator. Which of these labels is correct depends in many cases on the context, such as how a chord resolves to the next one. While these differences are important from a musical perspective (they express a different experience of the same musical events), our model was trained to only predict pitch-class distributions. Therefore, in its current state, it cannot reproduce these distinctions, but we expect future versions to significantly improve in this respect. End of Table 1