Improving an Iterative Physics Solver Using a Direct Method

[Music] roblox is a platform for user-generated 3d worlds of the games and the physics engine is a core component of the platform's runtime you this presentation will be an overview of the mathematics and implementation of the physics server inside of the roblox engine this is the system that computes the motion of constraining the rigid bodies under the influence of forces since the creativity of roblox developers knows no bounds or physics engine is required to support complex sometimes very complex mechanism for example we wouldn't be able to simulate this scissor-lift using an existing physics library this requires a more robust type of a constraint solver and

this will be the topic of this talk there is another requirement for a physics library the support for distributed physics which we will not discuss here the physics engine solves an NCP a linear complementarity problem and it typically it uses an interactive solver such as the projective gauss seidel within a reasonable number of iterations the protecting gauss-seidel is not able to simulate the scissor lift leading to instabilities and explosions our main contribution is the LDL PGS solver it uses a direct method to invert the sub matrix corresponding to the equality constraints and helping the PGS this leads to much better behavior in the example of the scissor lift also it behaves much better when bodies with very different masses are interacting the high mass ratio problem the LPGA's solver has two phases the symbolic phase a numeric phase in the

symbolic phase which runs only once per mechanism its structure is analyzed and an optimized LDL decomposition program is generated and in the numeric phase which executes every physics frame the mechanism state is updated using both the PGS and the LDL decomposition program in the next section we will quickly cover the basics of a constraint solver a constraint solver deals with mechanical constraints such as a ball in a socket or a hinge it also sports other constraints including rods prismatic cylindrical angular limits positional limits and ropes and few other constraints that don't need to have mechanical equivalents let's define a constraint mathematically we have two or more rigid bodies represented here by green and blue ovals and their state is described by their

coordinates positions and orientations the orientations can be specified using either rotation matrices or quaternions then a constraint between these bodies is a differentiable function of the coordinates into a vector space of an arbitrary dimension D this vector space is called the constrained space and it's dimension is called a degree of the constraint we say that the bodies respect the constraint if the constraint function vanishes that their coordinates the function must have an additional property it must be regular on the zero locus that is the Jacobian has full rank or equivalently D is the number of degrees of freedom removed by the constraint how do we model a ball in a socket or a ball joint between two rigid bodies specify pivot locations in the local space of the two bodies the constraint functions then the difference between these two pivot position transform into world space constraint function vanishes when the

two are in the same location it seems that there is another way to specify this constrain by using the tip status function between the two people positions but only the first model is regular that is the jacobian has full rank this is how we model a hinge it's a ball in the socket together with an orthogonality relationship between a vector on one body and the plane on the other body the degree of this constrain is 5 now that we have a way of modeling constraints we are interested in the dynamics of the rigid bodies their motion therefore we need to derive the velocity relationships at the constraining poses this is obtained by taking the derivative of the constraint function this derivative can always be formulated as a product of some matrix J times the vector of velocities the matrix is called the Jacobian of C

as a result the jacobian is an operator that takes the body space velocities and produces the constraint space velocities since the time derivative of a function that's zero is zero and strained bodies produce no motion in the constrained space for example here's the ball in a socket constraint and its derivative which is the relative velocity of the two pivots and care has to be taken in taking the derivative of orientation matrices or quaternions these are expressed in terms of cross products with angular velocities but a cross product is a linear operation and can be expressed as a matrix multiplication as a result we can gather the velocities into a vector and the rest into the Jacobian matrix the Jacobian matrix has two parts corresponding to the two rigid bodies so far we have been looking at single constraints but in practice we want to model entire mechanism consisting of many constraints and bodies we can stack

all the constraints in the mechanism into a single function called the global constraint taking the derivative of the global constraint we get the global Jacobian matrix which maps the velocities of rigid bodies in the constraint space practice each constraint we only affect two or three rigid bodies which means that each row of the global Jacobian has only twelve or eighteen nonzero entries here is an example of a mechanism with four bodies and four constraints and this is how the Jacobian would look like each block of rows corresponds to a constraint and each block of columns corresponds to a body the most important thing to notice is that each row of the Jacobian contains twelve nonzero elements the sparse matrix we are now ready to describe the integration step we have the velocities at time zero and using the Euler integration step the velocities at time DT are given by this equation it

involves the inverse mass matrix which is a block diagonal matrix with inverse masses and inverse inertia matrices on the diagonal and the Lagrangian multipliers are so called the constraint impulses this last term comes from the Dalembert's principle it adds the internal forces of the constraints and this internal forces must be such that the resulting velocities produced no motion in the constrained space moving the terms around we get the following equation it's called the linearized constraint equation it involves the constraint matrix K which is the product of the Jacobian inverse mass matrix and the Jacobian transpose this matrix is symmetric and positive semi definite which means its eigenvalues are non-negative it is either invertible or it can be made invertible by adding any positive value to its diagonal this equation is valid only for equality

constraints other constraints like collisions or limits lead to a null CP which we will not cover here but the structure of the problem is very similar in the next section we will discuss methods to solve this equation additionally the constraint equation or the NCP is solved using the projected gauss seidel this algorithmics small corrections to constraints consecutively the impulses for collisions and friction and for inequality constraints are processed using the projected gauss seidel iterations where the Equality constraints use the Gauss Seidel iteration these iterations are repeated multiple times until the system converges usually games use a fixed number of iterations our approach changes the way the Equality constraints are processed inside the PGS instead of processing these one at a time we use a direct method to invert the entire sub

matrix corresponding to the Equality constraints this fits within the formalism of the block gauss seidel and the method we use to invert the matrix is called the sparse LDL decomposition but because of performance we are not able to correct the solution using the entire holonomic block constraint matrix every iteration instead we do it only in one iteration the last one or one before last a lot of iterations use normal gauss seidel on the equality constraints it turns out that this is enough to have a very good quality simulation and as a bonus it allows us to compute the LDL decomposition in parallel with most of the PGS iterations let's review the PGS in details it starts with an estimate of the solution which can be cached result of the previous frame or simply the zero vector and it now starts the iterations over the rows of the constraint matrix

first calculates the residual which is the measure of the error for the current row and then corrects the solution by adding the residual divided by the diagonal element of the constraint matrix on the current row optionally it then projects the newly calculated solution into the admissible domain it repeats this process for all rows multiple time the bottleneck of this algorithm is this dot product between the row of K and lambda it has an arbitrary number of nonzero terms there is a better way to calculate the product ki times lambda let's go back to the constraint equation and write it in terms of the Jacobian and the inverse mass matrix we can use the associativity of the matrix multiplication to isolate the Jacobian on the left side and the product of the inverse mass matrix a cobia transpose and lambda on the right side we call the right side the vector of velocity changes Delta B it's literally the

velocity changes of the rigid bodies then in the gauss seidel we can replace ki times lambda by ji times delta b this product has only twelve or eighteen nonzero terms for practical strains we have now an improved version of the gauss seidel for mechanical constraints the physics simulation community it's called the impulse solver but in mathematics it's the catch manage method the solution is initialized to an estimate and the vector of velocity changes is pre computed as before we're gonna iterate over all rows of the matrix correction to the solution is calculated with the residue are computed using the new method optionally the corrected solution is projected and the velocity changes are updated this new algorithm uses about 50 floating-point operations per row and a constant number of operations is good

for an algorithm we already talked about the block gauss seidel this is formalized using partitions of rows a partition is a grouping into subsets where each element is included in exactly one subset a partition of rows induces a block structure on the jacobian it's the original Jacobian but where the rows have been permitted and grouped together according to the partition the block structure on the Jacobian naturally induces a block structure of the constraint matrix note that the square blocks on the diagonal are symmetric matrices and the opposite of diagonal blocks are transposed of each other we will call this using N and E and stands for node matrix and E for edge matrix we will see later this comes from graph theory the Gauss Seidel or the impulse solver can be made to work with blocks set of individual rows the things to note

in the blocking pool solver we need to compute the inverse of a square matrix and I which might be expensive depending on the size of the block so the blocks over doesn't support projection so this can only be used with equality constraints in the blog version of the impulse over we need to evaluate the inverse of the non matrix times a vector or small dimensions this is as easy as inverting a matrix using naive method like hammers rule for example using the natural partition coming from grouping the rows back constraints the null matrices will have dimensions up to six by six and these are invertible because we assume the constraints are regular consider this idea what if we grouped all the equality constraints together in the partition but the subset PI 0 contains all the rows of the equality constraints and other subsets containing dividual rows and 0 is then a large

square matrix constraint matrix of equality constraints we denote this matrix using H H stands for holonomic we need now to evaluate the inverse of H times a vector but there are some obstacles to do that H may not be invertible this is solved by using regularization which adds compliance the mechanical constraints and makes the system invertible it may have large dimensions we deal with this by using sparse methods still it may have large tensor matrices where sparse methods won't help us we can increase parsecity by splitting bodies so how does regularization work H is symmetric positive semi definite which means the angin values or nuns negative this means that adding anything positive to the diagonal makes the matrix invertible of course adding absolute values is scale dependent and it's much better to add

epsilon times diagonal adding such a regularization term changes the mechanical system it makes the constraints compliant traditionally other PGS engines also use compliance to stabilize the solution using PGs on the regular ice system is equivalent to using PGS on the original system but damping the solutions after each raw update in bullets or all the e-library x' this term is called CFM or constraint force mixing let's now look at the sparsity of constraint matrices or a mechanism like this one which has a scissor lift on top of a platform with wheels the constraint matrix consists mostly of zeros this model has 50 constraints of the hinges the dimension of the matrix is 246 and it's density is 12% spatial SQL

DL processes this matrix in about 90,000 floating-point operations which on my machine takes about 70 microseconds so what is Olesky decomposition it says that I'm symmetric positive-definite matrix a and we decompose into a product l times D times L transpose where L is a lower triangular matrix with ones on the diagonal and these are diagonal matrix this is useful because both L inverse and D inverse can be efficiently evaluated as operators on a vector Oh for any vector B inverse times V can now be evaluated by applying the inverse of the decomposition factors computing the L inverse and its transpose times a vector and we're done using a method called back substitution which is very efficient and inverting D is trivial as

it involves only inverting scalars as suppose that a has a block structure coming from a partition of the roads and it has an LDL decomposition where D is a block diagonal and L is a lower triangular matrix with identity blocks on the diagonal and both these matrices inherit the block structure from a but why are we using a block LDL instead of simple scale version its performance block decomposition is much faster than in Scala decomposition as sparse block matrix requires less indexation than a scalar one next we will look into details of an algorithm that produces such decomposition there are at least two algorithms that produce LDL decompositions the Gaussian elimination and doolittle algorithm will use only Gaussian elimination in this presentation at the core of the blog Gaussian elimination there is an operation called

a shul complement then we have a symmetric block matrix where the blocks are denoted by n e e transpose and s we apply the following operation to this matrix we multiply the top row by e times n inverse and subtract from the bottom row as a result we get a matrix where the bottom left block has become 0 has been eliminated then we multiply the first column by n inverse times a transpose and subtract on the second column resulting the top right block being eliminated what we have done is to decompose the block matrix as a product of a lower triangular matrix block diagonal matrix and an upper triangular matrix which is the transpose of the lower triangular one and so we have obtained a block LDL decomposition the algorithm that produces this block LDL decomposition is the block Gaussian elimination which is

just a recursive computation of the short complement or our constraint matrix which is a symmetric matrix of a block structure will group blocks even further to simplify the notation we'll call this lower left column e star comma 0 and the bottom right block we'll call it s 0 we have this 2 by 2 block representation of the matrix and computing the show complement produces this decomposition we now apply recursively the shrewder complement here we have a matrix that has already been partially eliminated so we already have some blocks on the top left of the diagonal and underneath and to the right the elements have already been eliminated and this is the sub-matrix who will be working on it is the si minus 1 matrix from the previous iteration applying the short complement we get this decomposition where the matrices on the

left and right are called elementary matrices the algorithm underlying the shrew complement is divided in three steps number one invert the node matrix n number to compute the elementary matrix L and number three reduce the S matrix after the last iteration we will be left with a block diagonal matrix and the lower triangular matrix L that's the product of elementary matrices in the order that were produced and these are the factors of the LDL decomposition of our matrix we can now apply the inverse of H to a vector V and this is equivalent to applying the inverse of the factors of the LDL decomposition in the opposite order the inverse of the matrix L can be computed using its decomposition into the elementary matrices that we obtain during the Gaussian elimination the

inverse of each elementary matrix is simple to compute just negate the elements below the diagonal and don't evaluate the product of this matrices use this decomposition while applying to a vector and finally the inverse of D is the block diagonal matrix consisting of the inverses of the node matrices and eyes but these have already been computed during the Gaussian elimination let's turn our attention to the constraint matrix of equality constraints H the block structure comes from partitioning the rose by constraints so the diagonal blocks and eyes correspond to individual constraints and off diagonal blocks IJ's correspond to the interaction between two constraints if they have a rigid body in common usually a lot of the CI J's are 0 because not all constraints have a rigid bodies in common

and we use these zeros to help us reduce the number of operations during the Gaussian elimination yes absolutely this is possible and here is how we need to establish a relationship between space matrices and graphs translate the Gaussian elimination to a process on graphs and pack the sparse matrix data in memory in such a way that it optimizes the block Gaussian elimination let's dive into some graph theory first we have the body graph it's the graph where the nodes are rigid bodies and edges are constrained between them the dual of the body graph is the constraint graph it's also called the edge graph of the body graph because it's nodes are constraints and edges are common bodies between constraints as you can imagine the constraint graph is related to the constraint matrix what this means is that nodes correspond to diagonal blocks the node matrices and

edges correspond to the off diagonal blocks the H matrices and Gaussian elimination has an equivalent operational graph called graph elimination we'll dive into this later here's a simple example of a dune buggy with its body graph the nodes correspond to the following rigid bodies the main body the wheels the knuckles that joined the front wheels to the body in the steering mechanism and the steering rack taking the edge graph of the body graph we get the constraint graph this graph represents the structure of the constraint matrix of the equality constraints in this buggy let's turn our attention to the way the Gaussian elimination translates on to the constraint graph let's call the node and I corresponding to the node matrix the pivot the basic operation of taking the shul compliment at ni does the following to the graph one it adds edges between the neighbors of the people and two it eliminates all

edges with the pivot on this example graph let's select not 0 as the pivot and eliminate it this produces the graph where all neighbors of not 0 are joined by edges and edges with not 0 are removed note that this process creates four new edges meaning this fields in four zero blocks in the sparse matrix with nonzero values there is a better way if we permit our matrix in such a way that the rows corresponding to node one are at the top we can eliminate not one instead doing this produces this graph where no new edges have been created and therefore no zero elements of the sparse matrix have been filled in this is much better than eliminating nodes 0 so now on the resulting graph let's find another node to eliminate let's try not zero this creates a new edge but if we

eliminate node 3 there we know new edges created and the sparsity of the matrix will not increase continuing with this algorithm we produce an elimination order if it doesn't produce new edges it's called a perfect elimination order in this example we have a perfect elimination order and it's 1 3 0 to 4 but not all graphs have a perfect elimination order the objective of the elimination game is to find an order that minimizes the number of new edges created this is an np-complete problem but there are some very good heuristics here are two here is this the first one called minimum degree algorithm is a very popular method as it is very fast but it generates mediocre ordering it's mostly used for very large problems with millions of nodes

the second one minimum edge creation algorithm is more expensive but in general produces much better ordering this is our invention and is more appropriate for medium-sized problems where the graph does not exceed hundreds of nodes also we can afford running a more expensive algorithm because we only need to compute this once per mechanism and catch the results to reuse a refrain the minimum edge creation algorithm is very simple at its step it eliminates the people that will create the minimum number of new edges that is the people whose neighbors have the smallest number of missing connections between themselves it does require expensive service step but it is affordable on graphs of few hundreds of nodes by the end the graph elimination creates two lists an ordered sequence of nodes the people sequence and Forex people a

sequence of eliminated edges the edge sequences the elements of the edge sequences should be in the pivot order before proceeding let's enjoy this animation of the elimination of the scissor lift constrain graph with a graph elimination enhance let's go back to our sparse matrix this is symmetric we only store the blocks on and below the diagonal the node matrices on the diagonal are in the people sequence order and the H matrices appear below the diagonal also in the edge sequence order the light blue blocks are those that were not present in the original matrix but were filled in by the elimination these are initially zero but will acquire nonzero values during the process of decomposition

the colored blocks represent all blocks that we'll need for Gaussian elimination and we can pack this contiguously in memory we started blocks in column major order in memory but we store the elements of each block in row major order this memory layout is very important for efficient processing now let's look in details at the step of Gaussian elimination on the sparse matrix we have the Peapod matrix we LDL decomposing place using dense LDL which is always faster than using the Kramer's rule we compute the L matrices by multiplying the H matrices below the pivot matrix in0 inverse we start the L matrices in the temporary buffer this process can be done very efficiently the way the edge matrices are packed in memory it can be interpreted as one block matrix and

multiply by the inverse of the pivot we don't need to keep track of the block structure next we reduce this is the most expensive part of the process as it costs proportionally to the square of the height of the edge column again we don't need to keep track of the block structure above E&L matrices but we do need to have an index section of the target location in s not then finally we copy the L matrices from the temporary buffer over the age matrices that we will no longer use and we repeat this process with anyone at the or decompose n1 compute L matrices reduce and replace a with L matrices after repeating this process on all the pivots we are left with the block LDL decomposition with the L matrices replacing the edge matrices and the D matrices replacing the node matrices and

these node matrices are already LDL decomposed so it's easy to compute their inverse operator we have one key ingredient still missing indexation of the reduction step the results of the reduction step which consists of multiplying the edge matrices with the L matrices are being scattered over the rest of the sparse matrix the scattering needs to be indexed this indexation is pre computed four pairs of indices in two rows of E and L and it is called the reduction scattering indexation RSI it can use either pointers or indices of both node and edge matrices live in the same workspace this array can become quite large and using indices rather than pointers is better because it requires less memory for two bytes versus eight bytes for pointers on 64-bit CPUs

why is the performance of the block LDL better than the regular LDL let D be the width of an edge column which is the degree of the corresponding constraint or the pivot and H be its height we assume that H is bigger than D which is true most of the time for complex mechanism then we estimate the number of operations for the three phases of the show compliment and we estimate the memory access these are the dominant terms now modern processors are good at floating-point operations but are bad at memory access and if we compute the operation count over the memory access we get D as a consequence the best most efficient results are obtained with larger D as long as it stays smaller than H you

let's look at the implementation of the reduction step which represents about 70 percent of the processing time it's fairly simple we loop over the raw indices of both matrices and compute the dot product the result is subtracted from the locations index by the ASI under the degree parameter which is the degree of the pivot constraint is passed as a template argument this allows unrolling the dot product underlying the matrix multiplication which has very large impact on performance you now here is the implementation of a single pivot elimination we have the three steps here the composition of the node matrix which takes about five percent of the time the computation of L which takes about twenty five percent of time and should also be on roll over the degree and the reduction step

and here is how you invoke this function by using a sweet statement on the degree or dimension of the pivot the switch should handle all degrees I'll do those up to six are the most common and should be allowed to enroll their internal computations for performance the last topic we need to cover is about large dense sub matrices where space methods can't help let's look at this example of the track bagel this model has three components the main body and the two tracks the track consists of 40 plates connected with hinge constraint the constraint might risk as dimension 200 with 7.5% density it's nonzero elements are mostly very close to the diagonal such a system is easily decomposed using LDL sparse LDL processes this matrix in 43,000 operations in fact it has linear

complexity in the number of constraints in the track on the other hand the body of the vehicle has a totally different constraint matrix it has only 20 constraints which are mostly hinges and it's dimension is 88 but the density of the constraint matrix is 100 percent it takes 250,000 operations to decompose the complexity here is n cube in the number of wheels this is a problem which we have a solution called body shattering you all the shuttering we transform the system into a system of higher dimension 142 but with much lower than density of 17% here DL can decompose this in only 35,000 operations it turns out that this number of operations is linear in the number of wheels we went from n cube complexity to linear complexity and this is how the constraint matrix looks like after shattering

so then sub-matrices appear where we have a rigid body with many constraints like the body of the tank here we introduce a method called body shattering it finds bodies with many constraints where the sum of the degrees should be more than 20 splits the rigid body into smaller shards joins the shards using rigid joints and equally distribute the constraint over this chart you suppose this is the body graph of a mechanism with a central body overloaded with constraints and this is the corresponding constraint graph you can see that it is a complete graph The Associated sparse matrix is dense we get this rigid body graph by shattering the central body into three pieces and joining them using rigid joints the corresponding constraint graph has more loads but less edges and constraining matrix is sparser are more efficient to

LDL decompose since we're now near the conclusion of this talk let's put all the ingredients of the LDL PG s solver together for every mechanism we have two phases the symbolic phase that needs to be executed only once and in numeric phase that runs every phase extreme update the symbolic phase takes us an input the structure of the mechanism which is the connectivity data between bodies and constraints and produces the elimination sequence or the LDL program the numeric face takes the current state of the mechanism and executes the PGS solver in parallel with the LDL decomposition of the constraint matrix then after joining the two threads executes one step of the block gauss seidel using the LDL decomposed matrix to correct the constraint impulses there is a profile capture of a solver update

we can see the PGS solver running here and in parallel the LDL decomposition of the holonomic constraints finally the decomposed matrix is applied to the impulses here this is it I hope this talk will help you build a great physics engine thanks for watching you you