SDPSolver.cpp
// Compute the quantities needed to solve the Schur complement
// equation
//
// {{S, -B}, {B^T, 0}} . {dx, dy} = {r, s}
//
// (where S = SchurComplement, B = FreeVarMatrix), using the method
// described in the manual:
//
// - Compute S using BilinearPairingsXInv and BilinearPairingsY.
//
// - Stabilize S by adding a low-rank update S' = S + U U^T and
//   compute the Cholesky decomposition S' = L' L'^T.
//
// - Form B' = (B U) and compute
//
//   - SchurOffDiagonal = L'^{-1} B
//   - L'^{-1} U (this is stored implicitly in the stabilizeBlock* variables)
//   - Q = (L'^{-1} B')^T (L'^{-1} B') - {{0, 0}, {0, 1}}
//
// - Compute the LU decomposition of Q.
//
// This data is sufficient to efficiently solve the above equation for
// a given r,s.
//
// Inputs:
// - BilinearPairingsXInv, BilinearPairingsY (these are members of
//   SDPSolver, but we include them as arguments to emphasize that
//   they must be computed first)
// - choleskyStabilizeThreshold: the real constant \theta used in
//   stabilizing S
// Workspace (members of SDPSolver which are modified by this method
// and not used later):
// - SchurComplement
// - schurStabilizeIndices
// - schurStabilizeLambdas
// Outputs (members of SDPSolver which are modified by this method and
// used later):
// - SchurComplementCholesky
// - SchurOffDiagonal
// - stabilizeBlockIndices
// - stabilizeBlockUpdateRow
// - stabilizeBlockUpdateColumn
// - stabilizeBlocks
// - Q, Qpivots
//
Compute the quantities needed to solve the Schur complement
equation



(where S = SchurComplement, B = FreeVarMatrix), using the method
described in the manual.



Stabilize S by adding a low-rank update S' = S + U U^T and
compute the Cholesky decomposition S' = L' L'^T.

Form B' = (B U) and compute



Compute the LU decomposition of Q.

This data is sufficient to efficiently solve the above equation for
a given r,s.

BilinearPairingsXInv, BilinearPairingsY (these are members of
SDPSolver, but we include them as arguments to emphasize that
they must be computed first)
choleskyStabilizeThreshold: the real constant \theta used in
stabilizing S
Workspace (members of SDPSolver which are modified by this method
and not used later):
SchurComplement
schurStabilizeIndices
schurStabilizeLambdas
Outputs (members of SDPSolver which are modified by this method and
used later):
SchurComplementCholesky
SchurOffDiagonal
stabilizeBlockIndices
stabilizeBlockUpdateRow
stabilizeBlockUpdateColumn
stabilizeBlocks
Q, Qpivots
Compute the quantities needed to solve the Schur complement equation ( where S = SchurComplement B = FreeVarMatrix ) using the method described in the manual . Stabilize S by adding a low - rank update S ' = S + U U^T and compute the Cholesky decomposition S ' = L ' L ' ^T . Form B ' = ( B U ) and compute Compute the LU decomposition of Q . This data is sufficient to efficiently solve the above equation for a given r s . BilinearPairingsXInv BilinearPairingsY ( these are members of SDPSolver but we include them as arguments to emphasize that they must be computed first ) choleskyStabilizeThreshold : the real constant \ theta used in stabilizing S Workspace ( members of SDPSolver which are modified by this method and not used later ) : SchurComplement schurStabilizeIndices schurStabilizeLambdas Outputs ( members of SDPSolver which are modified by this method and used later ) : SchurComplementCholesky SchurOffDiagonal stabilizeBlockIndices stabilizeBlockUpdateRow stabilizeBlockUpdateColumn stabilizeBlocks Q Qpivots
void SDPSolver::initializeSchurComplementSolver(const BlockDiagonalMatrix &BilinearPairingsXInv,
                                                const BlockDiagonalMatrix &BilinearPairingsY,
                                                const Real &choleskyStabilizeThreshold) {
  computeSchurComplement(sdp, BilinearPairingsXInv, BilinearPairingsY, SchurComplement);
     L' L'^T = S' = SchurComplement + U U^T
   Here, the 'update' matrix U has columns given by
     U = ( Lambda_{p_1} e_{p_1}, ..., Lambda_{p_M} e_{p_M} )
   where e_p is a unit vector in the p-th direction and the
   Lambda_{p_m} are constants. The p_i are given by
   schurStabilizeIndices and the corresponding Lambda_i by
   schurStabilizeLambdas.
  choleskyDecompositionStabilized(SchurComplement, SchurComplementCholesky,
                                  schurStabilizeIndices,
                                  schurStabilizeLambdas,
                                  cast2double(choleskyStabilizeThreshold));
   SchurOffDiagonal = L'^{-1} FreeVarMatrix
  SchurOffDiagonal.copyFrom(sdp.FreeVarMatrix);
  blockMatrixLowerTriangularSolve(SchurComplementCholesky, SchurOffDiagonal);
   Next we compute L'^{-1} U, which is stored implicitly in terms of
   its nonzero submatrices in the stabilizeBlock* variables.
   total number of columns in the off-diagonal part B' = (B U)
   (currently just B; will accumulate the rest shortly)
  int offDiagonalColumns = SchurOffDiagonal.cols;
  stabilizeBlockIndices.resize(0);
  stabilizeBlockUpdateRow.resize(0);
  stabilizeBlockUpdateColumn.resize(0);
   j runs over blocks of SchurComplement
  for (unsigned int j = 0; j < SchurComplement.blocks.size(); j++) {
    if (schurStabilizeIndices[j].size() > 0) {
       the j-th block of S contains stabilized directions. We have a
       block submatrix
         U_j = stabilizeBlocks[j]
       of U for each such j.
      stabilizeBlockIndices.push_back(j);
       index of the first stabilized direction within the j-th block
      int startIndex = schurStabilizeIndices[j][0];
       set the dimensions of U_j 
      stabilizeBlocks[j].resize(SchurComplement.blocks[j].rows - startIndex,
                                schurStabilizeIndices[j].size());
       set U_j = 0
      stabilizeBlocks[j].setZero();
       for each column of U_j add Lambda_p in the row (p - startIndex)
      for (unsigned int c = 0; c < schurStabilizeIndices[j].size(); c++) {
        int r = schurStabilizeIndices[j][c] - startIndex;
        stabilizeBlocks[j].elt(r, c) = schurStabilizeLambdas[j][c];
      }
       append the row of U corresponding to the top-left of U_j
      stabilizeBlockUpdateRow.push_back(SchurComplement.blockStartIndices[j] + startIndex);
       append the column of U corresponding to the top-left of U_j
      stabilizeBlockUpdateColumn.push_back(offDiagonalColumns);
       update the number of off-diagonal columns
      offDiagonalColumns += stabilizeBlocks[j].cols;
    }
  }
   Set U = L'^{-1} U
   We do this by modifying the blocks U_j = stabilizeBlocks[j]
   in-place, multiplying by the inverse of the appropriate submatrix
   of SchurComplementCholesky.  We henceforth refer to L'^{-1} U as
   V to avoid confusion.
  #pragma omp parallel for schedule(dynamic)
  for (unsigned int j = 0; j < stabilizeBlockIndices.size(); j++) {
    int b = stabilizeBlockIndices[j];
    int startIndex = schurStabilizeIndices[b][0];
    Rtrsm("Left", "Lower", "NoTranspose", "NonUnitDiagonal",
          stabilizeBlocks[b].rows, stabilizeBlocks[b].cols, 1,
          &SchurComplementCholesky.blocks[b].elt(startIndex, startIndex),
          SchurComplementCholesky.blocks[b].rows,
          &stabilizeBlocks[b].elt(0, 0),
          stabilizeBlocks[b].rows);
  }
   Next, we compute
     Q = (L'^{-1} B')^T (L'^{-1} B') - {{0, 0}, {0, 1}}
   Where B' = (B U).  We think of Q as containing four blocks called
   Upper/Lower-Left/Right.
   Set the dimensions of Q
  Q.resize(offDiagonalColumns, offDiagonalColumns);
  Q.setZero();
   Here, SchurOffDiagonal = L'^{-1} B.
   UpperLeft(Q) = SchurOffDiagonal^T SchurOffDiagonal
  matrixSquareIntoBlock(SchurOffDiagonal, Q, 0, 0);
   Here, stabilizeBlocks contains the blocks of V = L'^{-1} U.
   LowerRight(Q) = V^T V - 1
  for (unsigned int j = 0; j < stabilizeBlockIndices.size(); j++) {
    int b = stabilizeBlockIndices[j];
    int c = stabilizeBlockUpdateColumn[j];
    matrixSquareIntoBlock(stabilizeBlocks[b], Q, c, c);
     subtract the identity matrix from this block
    for (int i = c; i < c + stabilizeBlocks[b].cols; i++)
      Q.elt(i, i) -= 1;
  }
   LowerLeft(Q) = V^T SchurOffDiagonal
  # pragma omp parallel for schedule(dynamic)
  for (unsigned int j = 0; j < stabilizeBlockIndices.size(); j++) {
    int b = stabilizeBlockIndices[j];
    int p = stabilizeBlockUpdateRow[j];
    int r = stabilizeBlockUpdateColumn[j];
    Rgemm("Transpose", "NoTranspose",
          stabilizeBlocks[b].cols,
          SchurOffDiagonal.cols,
          stabilizeBlocks[b].rows,
          1,
          &stabilizeBlocks[b].elements[0],
          stabilizeBlocks[b].rows,
          &SchurOffDiagonal.elt(p, 0),
          SchurOffDiagonal.rows,
          0,
          &Q.elt(r, 0),
          Q.rows);
  }
   UpperRight(Q) = LowerLeft(Q)^T
  # pragma omp parallel for schedule(static)
  for (int c = 0; c < SchurOffDiagonal.cols; c++)
    for (int r = SchurOffDiagonal.cols; r < Q.rows; r++)
      Q.elt(c, r) = Q.elt(r, c);
   Resize Qpivots appropriately and compute the LU decomposition of Q
  Qpivots.resize(Q.rows);
  LUDecomposition(Q, Qpivots);
}
void SDPSolver :: initializeSchurComplementSolver ( const BlockDiagonalMatrix & BilinearPairingsXInv , const BlockDiagonalMatrix & BilinearPairingsY , const Real & choleskyStabilizeThreshold ) { computeSchurComplement ( sdp , BilinearPairingsXInv , BilinearPairingsY , SchurComplement ) ; choleskyDecompositionStabilized ( SchurComplement , SchurComplementCholesky , schurStabilizeIndices , schurStabilizeLambdas , cast2double ( choleskyStabilizeThreshold ) ) ; SchurOffDiagonal . copyFrom ( sdp . FreeVarMatrix ) ; blockMatrixLowerTriangularSolve ( SchurComplementCholesky , SchurOffDiagonal ) ; int offDiagonalColumns = SchurOffDiagonal . cols ; stabilizeBlockIndices . resize ( 0 ) ; stabilizeBlockUpdateRow . resize ( 0 ) ; stabilizeBlockUpdateColumn . resize ( 0 ) ; for ( unsigned int j = 0 ; j < SchurComplement . blocks . size ( ) ; j ++ ) { if ( schurStabilizeIndices [ j ] . size ( ) > 0 ) { stabilizeBlockIndices . push_back ( j ) ; int startIndex = schurStabilizeIndices [ j ] [ 0 ] ; stabilizeBlocks [ j ] . resize ( SchurComplement . blocks [ j ] . rows - startIndex , schurStabilizeIndices [ j ] . size ( ) ) ; stabilizeBlocks [ j ] . setZero ( ) ; for ( unsigned int c = 0 ; c < schurStabilizeIndices [ j ] . size ( ) ; c ++ ) { int r = schurStabilizeIndices [ j ] [ c ] - startIndex ; stabilizeBlocks [ j ] . elt ( r , c ) = schurStabilizeLambdas [ j ] [ c ] ; } stabilizeBlockUpdateRow . push_back ( SchurComplement . blockStartIndices [ j ] + startIndex ) ; stabilizeBlockUpdateColumn . push_back ( offDiagonalColumns ) ; offDiagonalColumns += stabilizeBlocks [ j ] . cols ; } } #pragma  omp parallel for schedule(dynamic) for ( unsigned int j = 0 ; j < stabilizeBlockIndices . size ( ) ; j ++ ) { int b = stabilizeBlockIndices [ j ] ; int startIndex = schurStabilizeIndices [ b ] [ 0 ] ; Rtrsm ( " " , " " , " " , " " , stabilizeBlocks [ b ] . rows , stabilizeBlocks [ b ] . cols , 1 , & SchurComplementCholesky . blocks [ b ] . elt ( startIndex , startIndex ) , SchurComplementCholesky . blocks [ b ] . rows , & stabilizeBlocks [ b ] . elt ( 0 , 0 ) , stabilizeBlocks [ b ] . rows ) ; } Q . resize ( offDiagonalColumns , offDiagonalColumns ) ; Q . setZero ( ) ; matrixSquareIntoBlock ( SchurOffDiagonal , Q , 0 , 0 ) ; for ( unsigned int j = 0 ; j < stabilizeBlockIndices . size ( ) ; j ++ ) { int b = stabilizeBlockIndices [ j ] ; int c = stabilizeBlockUpdateColumn [ j ] ; matrixSquareIntoBlock ( stabilizeBlocks [ b ] , Q , c , c ) ; for ( int i = c ; i < c + stabilizeBlocks [ b ] . cols ; i ++ ) Q . elt ( i , i ) -= 1 ; } # pragma  omp parallel for schedule(dynamic) for ( unsigned int j = 0 ; j < stabilizeBlockIndices . size ( ) ; j ++ ) { int b = stabilizeBlockIndices [ j ] ; int p = stabilizeBlockUpdateRow [ j ] ; int r = stabilizeBlockUpdateColumn [ j ] ; Rgemm ( " " , " " , stabilizeBlocks [ b ] . cols , SchurOffDiagonal . cols , stabilizeBlocks [ b ] . rows , 1 , & stabilizeBlocks [ b ] . elements [ 0 ] , stabilizeBlocks [ b ] . rows , & SchurOffDiagonal . elt ( p , 0 ) , SchurOffDiagonal . rows , 0 , & Q . elt ( r , 0 ) , Q . rows ) ; } # pragma  omp parallel for schedule(static) for ( int c = 0 ; c < SchurOffDiagonal . cols ; c ++ ) for ( int r = SchurOffDiagonal . cols ; r < Q . rows ; r ++ ) Q . elt ( c , r ) = Q . elt ( r , c ) ; Qpivots . resize ( Q . rows ) ; LUDecomposition ( Q , Qpivots ) ; }
Compute the quantities needed to solve the Schur complement
equation
Compute the quantities needed to solve the Schur complement equation
// compute SchurComplementCholesky = L', where // //   L' L'^T = S' = SchurComplement + U U^T // // Here, the 'update' matrix U has columns given by // //   U = ( Lambda_{p_1} e_{p_1}, ..., Lambda_{p_M} e_{p_M} ) // // where e_p is a unit vector in the p-th direction and the // Lambda_{p_m} are constants. The p_i are given by // schurStabilizeIndices and the corresponding Lambda_i by // schurStabilizeLambdas. //  // SchurOffDiagonal = L'^{-1} FreeVarMatrix // Next we compute L'^{-1} U, which is stored implicitly in terms of // its nonzero submatrices in the stabilizeBlock* variables. // total number of columns in the off-diagonal part B' = (B U) // (currently just B; will accumulate the rest shortly) // j runs over blocks of SchurComplement // the j-th block of S contains stabilized directions. We have a // block submatrix // //   U_j = stabilizeBlocks[j] // // of U for each such j. // index of the first stabilized direction within the j-th block // set the dimensions of U_j  // set U_j = 0 // for each column of U_j add Lambda_p in the row (p - startIndex) // append the row of U corresponding to the top-left of U_j // append the column of U corresponding to the top-left of U_j // update the number of off-diagonal columns // Set U = L'^{-1} U // // We do this by modifying the blocks U_j = stabilizeBlocks[j] // in-place, multiplying by the inverse of the appropriate submatrix // of SchurComplementCholesky.  We henceforth refer to L'^{-1} U as // V to avoid confusion. // Next, we compute // //   Q = (L'^{-1} B')^T (L'^{-1} B') - {{0, 0}, {0, 1}} // // Where B' = (B U).  We think of Q as containing four blocks called // Upper/Lower-Left/Right. // Set the dimensions of Q // Here, SchurOffDiagonal = L'^{-1} B. // // UpperLeft(Q) = SchurOffDiagonal^T SchurOffDiagonal // Here, stabilizeBlocks contains the blocks of V = L'^{-1} U. // // LowerRight(Q) = V^T V - 1 // subtract the identity matrix from this block // LowerLeft(Q) = V^T SchurOffDiagonal // UpperRight(Q) = LowerLeft(Q)^T // Resize Qpivots appropriately and compute the LU decomposition of Q