Skip to content

Commit

Permalink
Markowitz pivoting (#64)
Browse files Browse the repository at this point in the history 
* WIP

* more WIP

* WIP

* compilation

* more work: elimination with column permutations and Markowitz rule,
sort of working

* more tests

* cleanup

* improved pivot selection, bug fixes

* forward transformation on F

* f backward transformation

* the remaining forward and backward transformations

* more tests

* complete BTRAN and FTRAN

* extra unittest

* clone permutation matrices

* turn off logging

* allow one instance of the eliminator to be invoked multiple times

* WIP on replacing the old LUFactorization code to use the new method

* invV is properly computed

* proper inversion of the basis matrix

* refactoring to how invV is computed

* even more simplification

* compilation issues due to renaming

* renaming for clarity, couple of permutation indices issues

* more stable pivoting
bug fix
better dumping of lu factors

* allow external triggering of refresh basis

* oops

* bug fix

* some debugging infrastructure

* magic flag
  • Loading branch information
@guykatzz
guykatzz authored Jun 12, 2018
1 parent 1dd1191 commit 86bf67b
Show file tree
Hide file tree
Showing 27 changed files with 2,489 additions and 1,428 deletions.
1 change: 1 addition & 0 deletions src/basis_factorization/BasisFactorizationError.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@ class BasisFactorizationError : public Error
UNKNOWN_BASIS_FACTORIZATION_TYPE = 2,
CORRUPT_PERMUATION_MATRIX = 3,
CANT_INVERT_BASIS_BECAUSE_BASIS_ISNT_AVAILABLE = 4,
GAUSSIAN_ELIMINATION_FAILED = 5,
};

BasisFactorizationError( BasisFactorizationError::Code code ) :
Expand Down
36 changes: 18 additions & 18 deletions src/basis_factorization/ForrestTomlinFactorization.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -151,8 +151,8 @@ void ForrestTomlinFactorization::pushEtaMatrix( unsigned columnIndex, const doub

_workQ.resetToIdentity();
for ( unsigned i = indexOfChangedUColumn; i < _m - 1; ++i )
_workQ._ordering[i] = i + 1;
_workQ._ordering[_m - 1] = indexOfChangedUColumn;
_workQ._rowOrdering[i] = i + 1;
_workQ._rowOrdering[_m - 1] = indexOfChangedUColumn;
_workQ.invert( _invWorkQ );

/*
Expand All @@ -164,7 +164,7 @@ void ForrestTomlinFactorization::pushEtaMatrix( unsigned columnIndex, const doub
it is not upper triangular.
*/
for ( unsigned i = 0; i < _m; ++i )
_workVector[i] = column[_invQ._ordering[i]];
_workVector[i] = column[_invQ._rowOrdering[i]];

for ( unsigned i = 0; i < _m; ++i )
{
Expand All @@ -190,8 +190,8 @@ void ForrestTomlinFactorization::pushEtaMatrix( unsigned columnIndex, const doub
{
// Initialize to the non-modified bump value, according to the
// permutations. This is cell (m,i).
unsigned originalRow = _workQ._ordering[_m - 1];
unsigned originalCol = _workQ._ordering[i];
unsigned originalRow = _workQ._rowOrdering[_m - 1];
unsigned originalCol = _workQ._rowOrdering[i];
double bumpValue = _U[originalCol]->_column[originalRow];

for ( const auto &a : newAs )
Expand All @@ -204,7 +204,7 @@ void ForrestTomlinFactorization::pushEtaMatrix( unsigned columnIndex, const doub
*/

bumpValue += a->_value *
_U[_workQ._ordering[i]]->_column[_workQ._ordering[a->_column]];
_U[_workQ._rowOrdering[i]]->_column[_workQ._rowOrdering[a->_column]];
}

// If the bump value is zero, nothing needs to be done
Expand All @@ -224,7 +224,7 @@ void ForrestTomlinFactorization::pushEtaMatrix( unsigned columnIndex, const doub

if ( i < _m - 1 )
{
double diagonalValue = _U[_workQ._ordering[i]]->_column[_workQ._ordering[i]];
double diagonalValue = _U[_workQ._rowOrdering[i]]->_column[_workQ._rowOrdering[i]];
ASSERT( !FloatUtils::isZero( diagonalValue ) );
newA->_value = - bumpValue / diagonalValue;
}
Expand Down Expand Up @@ -274,26 +274,26 @@ void ForrestTomlinFactorization::pushEtaMatrix( unsigned columnIndex, const doub
// row-wise, so this means permuting invQ's rows.

for ( unsigned i = 0; i < _m; ++i )
_workOrdering[i] = _invWorkQ._ordering[_Q._ordering[i]];
_workOrdering[i] = _invWorkQ._rowOrdering[_Q._rowOrdering[i]];
for ( unsigned i = 0; i < _m; ++i )
_Q._ordering[i] = _workOrdering[i];
_Q._rowOrdering[i] = _workOrdering[i];

// Recompute invQ
_Q.invert( _invQ );

// Ai = _Q * Ai * _invQ, for all new A matrices
for ( const auto &a : newAs )
{
a->_row = _invQ._ordering[a->_row];
a->_column = _invQ._ordering[a->_column];
a->_row = _invQ._rowOrdering[a->_row];
a->_column = _invQ._rowOrdering[a->_column];
}

// Finally, append the new As to the list
_A.append( newAs );

// If the number of A matrices is too great, condense them.
if ( _A.size() > GlobalConfiguration::REFACTORIZATION_THRESHOLD )
refactorizeBasis();
obtainFreshBasis();
}

void ForrestTomlinFactorization::forwardTransformation( const double *y, double *x ) const
Expand Down Expand Up @@ -363,7 +363,7 @@ void ForrestTomlinFactorization::forwardTransformation( const double *y, double

// Multiply by inv(Q)
for ( unsigned i = 0; i < _m; ++i )
_workW[i] = _workVector[_invQ._ordering[i]];
_workW[i] = _workVector[_invQ._rowOrdering[i]];

/****
Step 2: Find x such that: Um....U1 * invQ * x = w
Expand All @@ -384,7 +384,7 @@ void ForrestTomlinFactorization::forwardTransformation( const double *y, double

// We are now left with invQ x = w (for our modified w). Multiply by Q and be done.
for ( unsigned i = 0; i < _m; ++i )
x[i] = _workW[_Q._ordering[i]];
x[i] = _workW[_Q._rowOrdering[i]];
}

void ForrestTomlinFactorization::backwardTransformation( const double *y, double *x ) const
Expand Down Expand Up @@ -416,7 +416,7 @@ void ForrestTomlinFactorization::backwardTransformation( const double *y, double
// Note: this is easier to do with a column-wise representation of Q,
// which is just the row-wise representation of invQ.
for ( unsigned i = 0; i < _m; ++i )
_workVector[i] = y[_invQ._ordering[i]];
_workVector[i] = y[_invQ._rowOrdering[i]];

// Eliminate the U's
for ( unsigned i = 0; i < _m; ++i )
Expand Down Expand Up @@ -446,7 +446,7 @@ void ForrestTomlinFactorization::backwardTransformation( const double *y, double
// Note: this is easier to do with a column-wise representation of invQ,
// which is just the row-wise representation of Q.
for ( unsigned i = 0; i < _m; ++i )
x[i] = _workVector[_Q._ordering[i]];
x[i] = _workVector[_Q._rowOrdering[i]];

// Mutiply by the As
for ( auto a = _A.rbegin(); a != _A.rend(); ++a )
Expand Down Expand Up @@ -683,7 +683,7 @@ void ForrestTomlinFactorization::makeExplicitBasisAvailable()
if ( explicitBasisAvailable() )
return;

refactorizeBasis();
obtainFreshBasis();
}

const double *ForrestTomlinFactorization::getBasis() const
Expand Down Expand Up @@ -859,7 +859,7 @@ void ForrestTomlinFactorization::dump() const
printf( "*** Done dumping FT factorization ***\n\n" );
}

void ForrestTomlinFactorization::refactorizeBasis()
void ForrestTomlinFactorization::obtainFreshBasis()
{
for ( unsigned column = 0; column < _m; ++column )
{
Expand Down
2 changes: 1 addition & 1 deletion src/basis_factorization/ForrestTomlinFactorization.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -91,7 +91,7 @@ class ForrestTomlinFactorization : public IBasisFactorization
/*
Obtain the basis matrix from the oracle and compute a fresh factorization
*/
void refactorizeBasis();
void obtainFreshBasis();

public:
/*
Expand Down
Loading

0 comments on commit 86bf67b

Please sign in to comment.