   LU decomposition: adding/removing rows/cols
  I'm playing around with physics, specifically statics (like stacking).  I've managed to express the contact forces from an n body collision (say, a pyramid of blocks) as a system of equations for n contact points.  I was then able to express it in terms of the familiar Ax = b.  Solving for x lets me find the contact forces necessary to resolve the stack.I managed to successfully implement LU decomposition using Crout's algorithm (with partial pivoting) as explained here.  I mostly understand what's going on and how it works, though I'm hardly comfortable with it.Decomposing the matrix is an O(n^3) operation.  I can then use the decomposed matrix to solve for x given any arbitrary b in O(n^2) time.  As long as the objects don't actually move, this means I can change the forces acting on the bodies and solve the new resulting contact forces in O(n^2) time.  So if, for instance, the bodies represent a bridge I can solve the resulting contact forces (in O(n^2) time) for when something like a train passes over the bridge, and its weight has to get distributed differently as it moves.Suppose I now want to add another row/column, to represent adding a block to the stack, or another strut to the bridge.  I already have the LU decomposed form from before I added the new body.  Is there a way to leverage that to construct the new LU decomposed matrix?  Especially if doing so can lower the O(n^3) cost of decomposing the new matrix to something more like O(n^2).Likewise, suppose I want to remove a row/column, representing a strut in the bridge breaking, for instance.   
  YES!  Let the new LU be {{L,0},{L',1}} and {{U,U'},{0,U0}}, where L and U are the old L and U.  Multiplying out gives {{LU, LU'},{L'U,L'U'+U0}}.  LU' corresponds to the new column (minus the last row element), L'U corresponds to the new row (minus last column), and L'U'+U0 the corner element.  You can leverage the old LU decomposition to find L' and U' in O(n^2) time (and in parallel if you want).  U0 comes once you have L'U'.As far as removing a row/column, the simplest method is to proceed with elimination as usual, but ignore the row you want to remove, and set the corresponding element to 0.   ;  Actually, to delete a row/column you can probably just delete the row/column from the LU decomposition.  (i haven;t checked this formally)  This is because when you multiply L*U, the row/columns you delete from the LU only contribute to the corresponding rows columns in the original matrix.  Also this maintains triangularity   ;  Thanks, though I'm not familiar with your notation.  What does {{L,0},{L',1}} mean? (You could also post latex if that makes it easier; I know how to parse it).   ;  What I meant is the (m+1)x(m+1) matrix  L   |   0  ----------   L'  |   1where L is your original L matrix (size m x m ), 0 represents a zero vector (size m x 1), L' an unknown update vector to the L matrix (size 1 x m), and 1 represents a scalar.  I was trying to show the augmented matrix representing the new lower triangular matrix, and likewise for the upper matrix.When I write {{LU, LU'},{L'U,L'U'+U0}}, I use the same (although terrible) notation, with the same sizes:   LU  (mxm) | LU' (mx1)   ----------   L'U (1xm) | L'U'+U0 (1x1)Note that this corresponds to your new matrix that you want to decompose (adding one row/column on the edges).  Say that your matrix looks like (where A is the original mxm)   A   |   g   ----------   h   |   j(same format).  then (thankfully), we already know that A = LU, by definition.  So LU' = g,L'U = h, L'U'+U0 = j.  You know g,h,j from your physics.  Then solve for U', L'; once you do this find U0.I hope this is clearer (I am new at the forum thing, so stuff is clearer in my head than on paper)   ;  Ah, ok I see now.  L' is a row vector, U' is a column vector, and you find the equations from the fact that LU = A before and after the augmentation.Are there any singular matrix type problems to look out for?  Like if U0 = 0, then the matrix is singular, I think.Likewise with numerical issues.  If you were to build up the initial LU decomposition one variable at a time like this, that would be the same thing as no pivoting, right?  So I'm guessing that as you add rows you're slowly going to become numerically unstable.  Is there a way to catch for that?   ;  There should not be singularity issues because the LU decomp is unique (once you have constrained the diagonal elements of L (or U).  So the LU does not exist only if A is singular anyway.  As for numerical issues, I haven't implemented this.  There is probably a way to do swapping with this method too, but I don't know what it is.   