diff options
Diffstat (limited to 'Eigen/src/SVD')
-rw-r--r-- | Eigen/src/SVD/UpperBidiagonalization.h | 307 |
1 files changed, 280 insertions, 27 deletions
diff --git a/Eigen/src/SVD/UpperBidiagonalization.h b/Eigen/src/SVD/UpperBidiagonalization.h index 587de37a5..2a55d550d 100644 --- a/Eigen/src/SVD/UpperBidiagonalization.h +++ b/Eigen/src/SVD/UpperBidiagonalization.h @@ -31,7 +31,7 @@ template<typename _MatrixType> class UpperBidiagonalization typedef typename MatrixType::Index Index; typedef Matrix<Scalar, 1, ColsAtCompileTime> RowVectorType; typedef Matrix<Scalar, RowsAtCompileTime, 1> ColVectorType; - typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0> BidiagonalType; + typedef BandMatrix<RealScalar, ColsAtCompileTime, ColsAtCompileTime, 1, 0, RowMajor> BidiagonalType; typedef Matrix<Scalar, ColsAtCompileTime, 1> DiagVectorType; typedef Matrix<Scalar, ColsAtCompileTimeMinusOne, 1> SuperDiagVectorType; typedef HouseholderSequence< @@ -61,6 +61,7 @@ template<typename _MatrixType> class UpperBidiagonalization } UpperBidiagonalization& compute(const MatrixType& matrix); + UpperBidiagonalization& computeUnblocked(const MatrixType& matrix); const MatrixType& householder() const { return m_householder; } const BidiagonalType& bidiagonal() const { return m_bidiagonal; } @@ -85,45 +86,297 @@ template<typename _MatrixType> class UpperBidiagonalization bool m_isInitialized; }; -template<typename _MatrixType> -UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix) + +/** \internal + * Helper routine for the block bidiagonal reduction. + * This function reduces to bidiagonal form the left \c rows x \a blockSize vertical + * and \a blockSize x \c cols horizontal panels of the \a A. The bottom-right block + * is left unchanged, and the Arices X and Y. + */ +template<typename MatrixType> +void upperbidiagonalization_inplace_helper(MatrixType& A, + typename MatrixType::RealScalar *diagonal, + typename MatrixType::RealScalar *upper_diagonal, + typename MatrixType::Index bs, + Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> > X, + Ref<Matrix<typename MatrixType::Scalar, Dynamic, Dynamic> > Y) { - Index rows = matrix.rows(); - Index cols = matrix.cols(); + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + typedef Ref<Matrix<Scalar, Dynamic, 1> > SubColumnType; + typedef Ref<Matrix<Scalar, 1, Dynamic>, 0, InnerStride<> > SubRowType; + typedef Ref<Matrix<Scalar, Dynamic, Dynamic> > SubMatType; - eigen_assert(rows >= cols && "UpperBidiagonalization is only for matrices satisfying rows>=cols."); + Index brows = A.rows(); + Index bcols = A.cols(); + + Scalar tau_u, tau_u_prev(0), tau_v; + + for(Index k = 0; k < bs; ++k) + { + Index remainingRows = brows - k; + Index remainingSizeInBlock = bs - k - 1; + Index remainingCols = bcols - k - 1; + + SubMatType X_k1( X.block(k,0, remainingRows,k) ); + SubMatType V_k1( A.block(k,0, remainingRows,k) ); + + // 1 - update the k-th column of A + SubColumnType v_k = A.col(k).tail(remainingRows); + v_k -= V_k1 * Y.row(k).head(k).adjoint(); + if(k) v_k -= X_k1 * A.col(k).head(k); + + // 2 - construct left Householder transform in-place + v_k.makeHouseholderInPlace(tau_v, diagonal[k]); + + if(k+1<bcols) + { + SubMatType Y_k ( Y.block(k+1,0, remainingCols, k+1) ); + SubMatType U_k1 ( A.block(0,k+1, k,remainingCols) ); + + // this eases the application of Householder transforAions + // A(k,k) will store tau_v later + A(k,k) = Scalar(1); + + // 3 - Compute y_k^T = tau_v * ( A^T*v_k - Y_k-1*V_k-1^T*v_k - U_k-1*X_k-1^T*v_k ) + { + SubColumnType y_k( Y.col(k).tail(remainingCols) ); + + // let's use the begining of column k of Y as a temporary vector + SubColumnType tmp( Y.col(k).head(k) ); + y_k.noalias() = A.block(k,k+1, remainingRows,remainingCols).adjoint() * v_k; // bottleneck + tmp.noalias() = V_k1.adjoint() * v_k; + y_k.noalias() -= Y_k.leftCols(k) * tmp; + tmp.noalias() = X_k1.adjoint() * v_k; + y_k.noalias() -= U_k1.adjoint() * tmp; + y_k *= numext::conj(tau_v); + } + + // 4 - update k-th row of A (it will become u_k) + SubRowType u_k( A.row(k).tail(remainingCols) ); + u_k = u_k.conjugate(); + { + u_k -= Y_k * A.row(k).head(k+1).adjoint(); + if(k) u_k -= U_k1.adjoint() * X.row(k).head(k).adjoint(); + } + + // 5 - construct right Householder transform in-placecols + u_k.makeHouseholderInPlace(tau_u, upper_diagonal[k]); + + // this eases the application of Householder transforAions + // A(k,k+1) will store tau_u later + A(k,k+1) = Scalar(1); + + // 6 - Compute x_k = tau_u * ( A*u_k - X_k-1*U_k-1^T*u_k - V_k*Y_k^T*u_k ) + { + SubColumnType x_k ( X.col(k).tail(remainingRows-1) ); + + // let's use the begining of column k of X as a temporary vectors + // note that tmp0 and tmp1 overlaps + SubColumnType tmp0 ( X.col(k).head(k) ), + tmp1 ( X.col(k).head(k+1) ); + + x_k.noalias() = A.block(k+1,k+1, remainingRows-1,remainingCols) * u_k.transpose(); // bottleneck + tmp0.noalias() = U_k1 * u_k.transpose(); + x_k.noalias() -= X_k1.bottomRows(remainingRows-1) * tmp0; + tmp1.noalias() = Y_k.adjoint() * u_k.transpose(); + x_k.noalias() -= A.block(k+1,0, remainingRows-1,k+1) * tmp1; + x_k *= numext::conj(tau_u); + tau_u = numext::conj(tau_u); + u_k = u_k.conjugate(); + } + + if(k>0) A.coeffRef(k-1,k) = tau_u_prev; + tau_u_prev = tau_u; + } + else + A.coeffRef(k-1,k) = tau_u_prev; + + A.coeffRef(k,k) = tau_v; + } - m_householder = matrix; + if(bs<bcols) + A.coeffRef(bs-1,bs) = tau_u_prev; - ColVectorType temp(rows); + // update A22 + if(bcols>bs && brows>bs) + { + SubMatType A22( A.bottomRightCorner(brows-bs,bcols-bs) ); + SubMatType A21( A.block(bs,0, brows-bs,bs) ); + SubMatType A12( A.block(0,bs, bs, bcols-bs) ); + Scalar tmp = A12(bs-1,0); + A12(bs-1,0) = 1; + A22.noalias() -= A21 * Y.topLeftCorner(bcols,bs).bottomRows(bcols-bs).adjoint(); + A22.noalias() -= X.topLeftCorner(brows,bs).bottomRows(brows-bs) * A12; + A12(bs-1,0) = tmp; + } +} + +/** \internal + * + * Implementation of a block-bidiagonal reduction. + * It is based on the following paper: + * The Design of a Parallel Dense Linear Algebra Software Library: Reduction to Hessenberg, Tridiagonal, and Bidiagonal Form. + * by Jaeyoung Choi, Jack J. Dongarra, David W. Walker. (1995) + * section 3.3 + */ +template<typename MatrixType, typename BidiagType> +void upperbidiagonalization_inplace_blocked(MatrixType& A, BidiagType& bidiagonal, + typename MatrixType::Index maxBlockSize=32, + typename MatrixType::Scalar* tempData = 0) +{ + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + typedef Block<MatrixType,Dynamic,Dynamic> BlockType; + + Index rows = A.rows(); + Index cols = A.cols(); + Index size = (std::min)(rows, cols); + + typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixType::MaxColsAtCompileTime,1> TempType; + TempType tempVector; + if(tempData==0) + { + tempVector.resize(cols); + tempData = tempVector.data(); + } + + Matrix<Scalar,Dynamic,Dynamic,ColMajor> X(rows,maxBlockSize), Y(cols, maxBlockSize); + X.setOnes(); + Y.setOnes(); + + Index blockSize = (std::min)(maxBlockSize,size); + + + + Index k = 0; + for (k = 0; k < size; k += blockSize) + { + Index bs = (std::min)(size-k,blockSize); // actual size of the block + Index tcols = cols - k - bs; // trailing columns + Index trows = rows - k - bs; // trailing rows + Index brows = rows - k; // rows of the block + Index bcols = cols - k; // columns of the block + + // partition the matrix A: + // + // | A00 A01 A02 | + // | | + // A = | A10 A11 A12 | + // | | + // | A20 A21 A22 | + // + // where A11 is a bs x bs diagonal block, + // and performs the bidiagonalization of A11, A21, A12, without updating A22. + // + // A22 will be updated in a second stage using matrices X and Y and level 3 operations: + // A22 -= V*Y^T - X*U^T + // where V and U contains the left and right Householder vectors + // + // Finally, the algorithm continue on the updated A22. + + + + // Let: + // | A11 A12 | + // B = | | + // | A21 A22 | + BlockType B = A.block(k,k,brows,bcols); + upperbidiagonalization_inplace_helper<BlockType>(B, + &(bidiagonal.template diagonal<0>().coeffRef(k)), + &(bidiagonal.template diagonal<1>().coeffRef(k)), + bs, + X.topLeftCorner(brows,bs), + Y.topLeftCorner(bcols,bs) + ); + } +} + +// Standard upper bidiagonalization without fancy optimizations +// This version should be faster for small matrix size +template<typename MatrixType, typename BidiagType> +void upperbidiagonalization_inplace_unblocked(MatrixType& mat, BidiagType& bidiagonal, + typename MatrixType::Scalar* tempData = 0) +{ + typedef typename MatrixType::Index Index; + typedef typename MatrixType::Scalar Scalar; + + Index rows = mat.rows(); + Index cols = mat.cols(); + Index size = (std::min)(rows, cols); + + typedef Matrix<Scalar,Dynamic,1,ColMajor,MatrixType::MaxRowsAtCompileTime,1> TempType; + TempType tempVector; + if(tempData==0) + { + tempVector.resize(rows); + tempData = tempVector.data(); + } for (Index k = 0; /* breaks at k==cols-1 below */ ; ++k) { Index remainingRows = rows - k; Index remainingCols = cols - k - 1; - // construct left householder transform in-place in m_householder - m_householder.col(k).tail(remainingRows) - .makeHouseholderInPlace(m_householder.coeffRef(k,k), - m_bidiagonal.template diagonal<0>().coeffRef(k)); - // apply householder transform to remaining part of m_householder on the left - m_householder.bottomRightCorner(remainingRows, remainingCols) - .applyHouseholderOnTheLeft(m_householder.col(k).tail(remainingRows-1), - m_householder.coeff(k,k), - temp.data()); + // construct left householder transform in-place in A + mat.col(k).tail(remainingRows) + .makeHouseholderInPlace(mat.coeffRef(k,k), bidiagonal.template diagonal<0>().coeffRef(k)); + // apply householder transform to remaining part of A on the left + mat.bottomRightCorner(remainingRows, remainingCols) + .applyHouseholderOnTheLeft(mat.col(k).tail(remainingRows-1), mat.coeff(k,k), tempData); if(k == cols-1) break; - - // construct right householder transform in-place in m_householder - m_householder.row(k).tail(remainingCols) - .makeHouseholderInPlace(m_householder.coeffRef(k,k+1), - m_bidiagonal.template diagonal<1>().coeffRef(k)); - // apply householder transform to remaining part of m_householder on the left - m_householder.bottomRightCorner(remainingRows-1, remainingCols) - .applyHouseholderOnTheRight(m_householder.row(k).tail(remainingCols-1).transpose(), - m_householder.coeff(k,k+1), - temp.data()); + + // construct right householder transform in-place in mat + mat.row(k).tail(remainingCols) + .makeHouseholderInPlace(mat.coeffRef(k,k+1), bidiagonal.template diagonal<1>().coeffRef(k)); + // apply householder transform to remaining part of mat on the left + mat.bottomRightCorner(remainingRows-1, remainingCols) + .applyHouseholderOnTheRight(mat.row(k).tail(remainingCols-1).transpose(), mat.coeff(k,k+1), tempData); } +} + +template<typename _MatrixType> +UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::computeUnblocked(const _MatrixType& matrix) +{ + Index rows = matrix.rows(); + Index cols = matrix.cols(); + + eigen_assert(rows >= cols && "UpperBidiagonalization is only for Arices satisfying rows>=cols."); + + m_householder = matrix; + + ColVectorType temp(rows); + + upperbidiagonalization_inplace_unblocked(m_householder, m_bidiagonal, temp.data()); + + MatrixType A = matrix; + BidiagonalType B(cols,cols); + upperbidiagonalization_inplace_blocked(A, B, 8); + + std::cout << (m_householder-A)/*.maxCoeff()*/ << "\n\n"; + std::cout << m_householder << "\n\n" + << m_bidiagonal.template diagonal<0>() << "\n" + << m_bidiagonal.template diagonal<1>() << "\n\n"; + std::cout << A << "\n\n" + << B.template diagonal<0>() << "\n" + << B.template diagonal<1>() << "\n\n"; + + m_isInitialized = true; + return *this; +} + +template<typename _MatrixType> +UpperBidiagonalization<_MatrixType>& UpperBidiagonalization<_MatrixType>::compute(const _MatrixType& matrix) +{ + Index rows = matrix.rows(); + Index cols = matrix.cols(); + + eigen_assert(rows >= cols && "UpperBidiagonalization is only for Arices satisfying rows>=cols."); + + m_householder = matrix; + upperbidiagonalization_inplace_blocked(m_householder, m_bidiagonal); + m_isInitialized = true; return *this; } |