From 47cb4fb2ed8feff34d1cefaa7645861b63f91fc5 Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Tue, 2 May 2017 20:17:19 +0200
Subject: [PATCH 01/10] Implement gemv

---
 src/internal_utils.rs | 43 +++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 41 insertions(+), 2 deletions(-)

diff --git a/src/internal_utils.rs b/src/internal_utils.rs
index dd2818f..5a4d8f5 100644
--- a/src/internal_utils.rs
+++ b/src/internal_utils.rs
@@ -1,6 +1,6 @@
 use matrix::{BaseMatrix, BaseMatrixMut};
 use libnum::{Zero, Num};
-use utils::in_place_vec_bin_op;
+use utils::{dot, in_place_vec_bin_op};
 
 pub fn nullify_lower_triangular_part<T, M>(matrix: &mut M)
     where T: Zero, M: BaseMatrixMut<T> {
@@ -20,6 +20,30 @@ pub fn nullify_upper_triangular_part<T, M>(matrix: &mut M)
     }
 }
 
+/// Given a vector `x` and a `m x n` matrix `A`, compute
+/// `y = A x`.
+///
+/// This is a stopgap solution until we have a more proper
+/// BLIS/BLAS-like API.
+pub fn gemv<T, M>(a: &M, x: &[T], y: &mut [T])
+    where M: BaseMatrix<T>, T: Num + Copy
+{
+    let m = a.rows();
+    let n = a.cols();
+
+    assert!(x.len() == n, "A and x must be dimensionally compatible.");
+    assert!(y.len() == m, "A and y must be dimensionally compatible.");
+
+    for element in y.iter_mut() {
+        *element = T::zero();
+    }
+
+    for i in 0 .. m {
+        let a_i = a.row(i).raw_slice();
+        y[i] = dot(a_i, x);
+    }
+}
+
 /// Given a vector `x` and a `m x n` matrix `A`, compute
 /// `y = A^T x`.
 ///
@@ -88,7 +112,7 @@ mod tests {
 
     use super::nullify_lower_triangular_part;
     use super::nullify_upper_triangular_part;
-    use super::transpose_gemv;
+    use super::{gemv, transpose_gemv};
     use super::ger;
 
     #[test]
@@ -117,6 +141,21 @@ mod tests {
         ]);
     }
 
+    #[test]
+    fn gemv_examples() {
+        {
+            let a = matrix![3.0, 4.0;
+                            5.0, 2.0;
+                            3.0, 1.0];
+            let x = vec![2.0, 3.0];
+            let mut y = vec![0.0; 3];
+            gemv(&a, &x, &mut y);
+
+            let y = Vector::new(y);
+            assert_vector_eq!(y, vector![18.0, 16.0, 9.0]);
+        }
+    }
+
     #[test]
     fn transpose_gemv_examples() {
         {

From f544a32c8fa601384dcdd1ec4e317d7ba6b32291 Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Tue, 2 May 2017 20:21:51 +0200
Subject: [PATCH 02/10] Implement HouseholderReflection right multiply

---
 src/matrix/decomposition/householder.rs | 54 ++++++++++++++++++++++++-
 1 file changed, 53 insertions(+), 1 deletion(-)

diff --git a/src/matrix/decomposition/householder.rs b/src/matrix/decomposition/householder.rs
index 4b26b20..6916df9 100644
--- a/src/matrix/decomposition/householder.rs
+++ b/src/matrix/decomposition/householder.rs
@@ -1,6 +1,7 @@
 use matrix::{Matrix, BaseMatrix, BaseMatrixMut, Column, ColumnMut};
 use vector::Vector;
 use utils;
+use internal_utils::{gemv, transpose_gemv, ger};
 
 use libnum::Float;
 
@@ -100,7 +101,6 @@ impl<T: Float> HouseholderReflection<T> {
     pub fn buffered_left_multiply_into<M>(&self, matrix: &mut M, buffer: &mut [T])
         where M: BaseMatrixMut<T>
     {
-        use internal_utils::{transpose_gemv, ger};
         assert!(buffer.len() == matrix.cols());
 
         // Recall that the Householder reflection is represented by
@@ -129,6 +129,35 @@ impl<T: Float> HouseholderReflection<T> {
         ger(matrix, - self.tau, v, u);
     }
 
+    /// Right-multiplies the given matrix by this Householder reflection.
+    ///
+    /// More precisely, let `H` denote this Householder reflection matrix,
+    /// and let `A` be a dimensionally compatible matrix. Then
+    /// this function computes the product `AH` and stores the result
+    /// back in `A`.
+    ///
+    /// The user must provide a buffer of size `A.rows()` which is used
+    /// to store intermediate results.
+    pub fn buffered_right_multiply_into<M>(&self, matrix: &mut M, buffer: &mut [T])
+        where M: BaseMatrixMut<T>
+    {
+        assert!(buffer.len() == matrix.rows());
+
+        // See `buffered_left_multiply_into`. The implementation here
+        // is almost exactly the same, except we instead have:
+        // AH = A - τ (A v) vᵀ = A - τ u vᵀ,
+        // where u = Av.
+
+        let ref v = self.v.data();
+        let mut u = buffer;
+
+        // u = A v
+        gemv(matrix, v, u);
+
+        // A <- A - τ v uᵀ
+        ger(matrix, - self.tau, u, v);
+    }
+
     pub fn as_vector(&self) -> &Vector<T> {
         &self.v
     }
@@ -377,6 +406,29 @@ mod tests {
         assert_matrix_eq!(x, expected, comp = abs, tol = 1e-3);
     }
 
+    #[test]
+    fn householder_reflection_right_multiply() {
+        let mut x = matrix![ 0.0,  1.0,  2.0,  3.0;
+                             4.0,  5.0,  6.0,  7.0;
+                             8.0,  9.0, 10.0, 11.0;
+                            12.0, 13.0, 14.0, 15.0 ];
+
+        let h = HouseholderReflection {
+            tau: 1.0 / 15.0,
+            v: vector![1.0, 2.0, 3.0, 4.0]
+        };
+
+        let mut buffer = vec![0.0; 4];
+
+        h.buffered_right_multiply_into(&mut x, &mut buffer);
+
+        let expected = matrix![ -1.3333, -1.6667, - 2.0000, - 2.3333;
+                                 0.0000, -3.0000, - 6.0000, - 9.0000;
+                                 1.3333, -4.3333, -10.0000, -15.6667;
+                                 2.6667, -5.6667, -14.0000, -22.3333 ];
+        assert_matrix_eq!(x, expected, comp = abs, tol = 1e-3);
+    }
+
     #[test]
     fn householder_composition_left_multiply() {
         let storage = matrix![ 5.0,  3.0,  2.0;

From d0bfc316b4e1a24afe68dc7168a48f6961af6503 Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Fri, 5 May 2017 11:20:02 +0200
Subject: [PATCH 03/10] Implement is_upper_hessenberg

Necessary for testing the validity of the output
from Hessenberg decomposition.
---
 src/testsupport/constraints.rs | 77 +++++++++++++++++++++++++++++++++-
 src/testsupport/mod.rs         |  3 +-
 2 files changed, 77 insertions(+), 3 deletions(-)

diff --git a/src/testsupport/constraints.rs b/src/testsupport/constraints.rs
index 4078645..3ca58dd 100644
--- a/src/testsupport/constraints.rs
+++ b/src/testsupport/constraints.rs
@@ -30,10 +30,28 @@ pub fn is_upper_triangular<T, M>(m: &M) -> bool
                         .all(|x| x == &T::zero()))
 }
 
+/// Returns true if the matrix is upper Hessenberg, otherwise false.
+/// Note that if the matrix is not square, it cannot be
+/// upper Hessenberg, and so this function returns false.
+pub fn is_upper_hessenberg<T, M>(m: &M) -> bool
+    where T: Zero + PartialEq<T>, M: BaseMatrix<T>
+{
+    if m.rows() == m.cols() {
+        m.row_iter()
+         .enumerate()
+         .skip(1)
+         .all(|(i, row)| row.iter().take(i - 1)
+                                   .all(|x| x == &T::zero()))
+    } else {
+        false
+    }
+}
+
 #[cfg(test)]
 mod tests {
     use super::is_lower_triangular;
     use super::is_upper_triangular;
+    use super::is_upper_hessenberg;
     use matrix::Matrix;
 
     #[test]
@@ -122,18 +140,73 @@ mod tests {
         }
     }
 
+    #[test]
+    fn is_upper_hessenberg_valid_input() {
+        {
+            let x: Matrix<u32> = matrix![];
+            assert!(is_upper_hessenberg(&x));
+        }
+
+        {
+            let x = matrix![3];
+            assert!(is_upper_hessenberg(&x));
+        }
+
+        {
+            let x = matrix![3, 2;
+                            4, 0];
+            assert!(is_upper_hessenberg(&x));
+        }
+
+        {
+            let x = matrix![3, 2, 0;
+                            4, 0, 2;
+                            0, 1, 3];
+            assert!(is_upper_hessenberg(&x));
+        }
+
+        {
+            let x = matrix![3, 2, 0, 2;
+                            4, 0, 2, 1;
+                            0, 1, 3, 4;
+                            0, 0, 2, 1];
+            assert!(is_upper_hessenberg(&x));
+        }
+    }
+
+    #[test]
+    fn is_upper_hessenberg_invalid_input() {
+        {
+            let x = matrix![3, 2, 0;
+                            4, 0, 2;
+                            3, 1, 3];
+            assert!(!is_upper_hessenberg(&x));
+        }
+
+        {
+            let x = matrix![3, 2, 0, 5;
+                            4, 0, 2, 2;
+                            3, 1, 3, 1;
+                            0, 3, 2, 0];
+            assert!(!is_upper_hessenberg(&x));
+        }
+    }
+
     quickcheck! {
         fn property_zero_is_lower_triangular(m: usize, n: usize) -> bool {
             let x = Matrix::<u32>::zeros(m, n);
             is_lower_triangular(&x)
         }
-    }
 
-    quickcheck! {
         fn property_zero_is_upper_triangular(m: usize, n: usize) -> bool {
             let x = Matrix::<u32>::zeros(m, n);
             is_upper_triangular(&x)
         }
+
+        fn property_zero_is_upper_hessenberg(n: usize) -> bool {
+            let x = Matrix::<u32>::zeros(n, n);
+            is_upper_hessenberg(&x)
+        }
     }
 }
 
diff --git a/src/testsupport/mod.rs b/src/testsupport/mod.rs
index 3277d1f..9356a0e 100644
--- a/src/testsupport/mod.rs
+++ b/src/testsupport/mod.rs
@@ -2,4 +2,5 @@
 //! tests involving data structures provided by rulinalg.
 mod constraints;
 
-pub use self::constraints::{is_lower_triangular, is_upper_triangular};
+pub use self::constraints::{is_lower_triangular, is_upper_triangular,
+                            is_upper_hessenberg};

From 832ff2e34ffeeafbdd725d41ed356b2467e9d5b0 Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Fri, 5 May 2017 13:47:13 +0200
Subject: [PATCH 04/10] Relax restrictions on buffer for Householder

We don't need the buffer to be a specific size,
it's sufficient that it's large enough to hold
the intermediate data. This simplifies usage
a little, since we don't need to truncate or
resize the buffer when it is reused for different
instances of HouseholderReflection.
---
 src/matrix/decomposition/householder.rs | 12 ++++++------
 src/matrix/decomposition/qr.rs          |  2 +-
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/src/matrix/decomposition/householder.rs b/src/matrix/decomposition/householder.rs
index 6916df9..7e3db4e 100644
--- a/src/matrix/decomposition/householder.rs
+++ b/src/matrix/decomposition/householder.rs
@@ -96,12 +96,12 @@ impl<T: Float> HouseholderReflection<T> {
     /// this function computes the product `HA` and stores the result
     /// back in `A`.
     ///
-    /// The user must provide a buffer of size `A.cols()` which is used
+    /// The user must provide a buffer of size at least `A.cols()` which is used
     /// to store intermediate results.
     pub fn buffered_left_multiply_into<M>(&self, matrix: &mut M, buffer: &mut [T])
         where M: BaseMatrixMut<T>
     {
-        assert!(buffer.len() == matrix.cols());
+        assert!(buffer.len() >= matrix.cols());
 
         // Recall that the Householder reflection is represented by
         // H = I - τ v vᵀ,
@@ -120,7 +120,7 @@ impl<T: Float> HouseholderReflection<T> {
         // Instead, we will use the provided buffer to hold the result of the
         // matrix-vector product.
         let ref v = self.v.data();
-        let mut u = buffer;
+        let mut u = &mut buffer[0 .. matrix.cols()];
 
         // u = A^T v
         transpose_gemv(matrix, v, u);
@@ -136,12 +136,12 @@ impl<T: Float> HouseholderReflection<T> {
     /// this function computes the product `AH` and stores the result
     /// back in `A`.
     ///
-    /// The user must provide a buffer of size `A.rows()` which is used
+    /// The user must provide a buffer of size at least `A.rows()` which is used
     /// to store intermediate results.
     pub fn buffered_right_multiply_into<M>(&self, matrix: &mut M, buffer: &mut [T])
         where M: BaseMatrixMut<T>
     {
-        assert!(buffer.len() == matrix.rows());
+        assert!(buffer.len() >= matrix.rows());
 
         // See `buffered_left_multiply_into`. The implementation here
         // is almost exactly the same, except we instead have:
@@ -149,7 +149,7 @@ impl<T: Float> HouseholderReflection<T> {
         // where u = Av.
 
         let ref v = self.v.data();
-        let mut u = buffer;
+        let mut u = &mut buffer[0..matrix.rows()];
 
         // u = A v
         gemv(matrix, v, u);
diff --git a/src/matrix/decomposition/qr.rs b/src/matrix/decomposition/qr.rs
index 31dd9f6..11261ab 100644
--- a/src/matrix/decomposition/qr.rs
+++ b/src/matrix/decomposition/qr.rs
@@ -206,7 +206,7 @@ impl<T> HouseholderQr<T> where T: Float {
             // gets shorter for each iteration, so we truncate the buffer
             // to the appropriate length.
             buffer.truncate(m - j);
-            multiply_buffer.truncate(bottom_right.cols());
+
             bottom_right.col(0).clone_into_slice(&mut buffer);
 
             let house = HouseholderReflection::compute(Vector::new(buffer));

From deaa23384816626e5d8d09df9692049764dfe6b8 Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Fri, 5 May 2017 17:34:39 +0200
Subject: [PATCH 05/10] Generalize HouseholderComposition with offset

The composition of Householder matrices for
QR decomposition and Upper Hessenberg decomposition
is very similar. We can exploit this to reuse virtually
all code we have written for HouseholderComposition in the
QR decomposition to implement the Hessenberg decomposition.
We do this by adding a "subdiagonal" parameter to
HouseholderComposition which controls the
subdiagonal to pull Householder reflectors from.
---
 src/matrix/decomposition/householder.rs | 38 +++++++++++++++----------
 src/matrix/decomposition/qr.rs          |  6 ++--
 2 files changed, 26 insertions(+), 18 deletions(-)

diff --git a/src/matrix/decomposition/householder.rs b/src/matrix/decomposition/householder.rs
index 7e3db4e..0458555 100644
--- a/src/matrix/decomposition/householder.rs
+++ b/src/matrix/decomposition/householder.rs
@@ -206,7 +206,11 @@ impl<T: Float> HouseholderReflection<T> {
 #[derive(Debug, Clone)]
 pub struct HouseholderComposition<'a, T> where T: 'a {
     storage: &'a Matrix<T>,
-    tau: &'a [T]
+    tau: &'a [T],
+
+    // The subdiagonal index determines the subdiagonal to pull Householder
+    // reflectors from. This is 0 for QR and 1 for Hessenberg.
+    subdiagonal: usize
 }
 
 /// Instantiates a HouseholderComposition with the given
@@ -217,12 +221,13 @@ pub struct HouseholderComposition<'a, T> where T: 'a {
 /// a HouseholderComposition by themselves, which is desirable
 /// because we want to have the freedom to change details
 /// of the internal representation if necessary.
-pub fn create_composition<'a, T>(storage: &'a Matrix<T>, tau: &'a [T])
+pub fn create_composition<'a, T>(storage: &'a Matrix<T>, tau: &'a [T], subdiagonal: usize)
     -> HouseholderComposition<'a, T>
 {
     HouseholderComposition {
             storage: storage,
-            tau: tau
+            tau: tau,
+            subdiagonal: subdiagonal
     }
 }
 
@@ -238,15 +243,17 @@ impl<'a, T> HouseholderComposition<'a, T> where T: Float {
         let n = self.storage.cols();
         let p = min(m, n);
         let q = matrix.cols();
+        let s = self.subdiagonal;
 
         assert!(matrix.rows() == m, "Matrix does not have compatible dimensions.");
 
         let mut house_buffer = Vec::with_capacity(m);
         let mut multiply_buffer = vec![T::zero(); q];
-        for j in (0 .. p).rev() {
-            house_buffer.resize(m - j, T::zero());
-            let storage_block = self.storage.sub_slice([j, j], m - j, n - j);
-            let mut matrix_block = matrix.sub_slice_mut([j, 0], m - j, q);
+        for j in (0 .. p.saturating_sub(s)).rev() {
+            house_buffer.resize(m - j - s, T::zero());
+
+            let storage_block = self.storage.sub_slice([j + s, j], m - j - s, n - j);
+            let mut matrix_block = matrix.sub_slice_mut([j + s, 0], m - j - s, q);
             let house = load_house_from_col(&storage_block.col(0),
                                             self.tau[j], house_buffer);
             house.buffered_left_multiply_into(&mut matrix_block,
@@ -266,6 +273,7 @@ impl<'a, T> HouseholderComposition<'a, T> where T: Float {
         let m = self.storage.rows();
         let n = self.storage.cols();
         let p = min(m, n);
+        let s = self.subdiagonal;
 
         assert!(k <= self.storage.rows(),
             "k cannot exceed m, the number of rows of Q");
@@ -281,12 +289,12 @@ impl<'a, T> HouseholderComposition<'a, T> where T: Float {
         // to reduce the number of operations
         // (note the size of the "q_k_block")
         let mut buffer = Vec::with_capacity(m);
-        let mut multiply_buffer = Vec::with_capacity(k);
-        for j in (0 .. min(p, k)).rev() {
-            buffer.resize(m - j, T::zero());
-            multiply_buffer.resize(k - j, T::zero());
-            let storage_block = self.storage.sub_slice([j, j], m - j, n - j);
-            let mut q_k_block = q_k.sub_slice_mut([j, j], m - j, k - j);
+        let mut multiply_buffer = vec![T::zero(); k];
+        for j in (0 .. min(p.saturating_sub(s), k)).rev() {
+            buffer.resize(m - j - s, T::zero());
+
+            let storage_block = self.storage.sub_slice([j + s, j], m - j - s, n - j);
+            let mut q_k_block = q_k.sub_slice_mut([j + s, j], m - j - s, k - j);
             let house = load_house_from_col(&storage_block.col(0),
                                             self.tau[j], buffer);
             house.buffered_left_multiply_into(&mut q_k_block,
@@ -443,7 +451,7 @@ mod tests {
         // let q = matrix![7.0/9.0, -28.0/45.0,   4.0/45.0;
         //                -4.0/9.0, - 4.0/ 9.0,   7.0/ 9.0;
         //                 4.0/9.0,  29.0/45.0,  28.0/45.0];
-        let composition = create_composition(&storage, &tau);
+        let composition = create_composition(&storage, &tau, 0);
 
         {
             // Square
@@ -489,7 +497,7 @@ mod tests {
                                2.0,  1.0,  3.0;
                               -2.0,  3.0, -2.0];
         let tau = vec![2.0/9.0, 1.0 / 5.0, 2.0];
-        let composition = create_composition(&storage, &tau);
+        let composition = create_composition(&storage, &tau, 0);
 
         // This corresponds to the following `Q` matrix
         let q = matrix![7.0/9.0, -28.0/45.0,   4.0/45.0;
diff --git a/src/matrix/decomposition/qr.rs b/src/matrix/decomposition/qr.rs
index 11261ab..412c1d6 100644
--- a/src/matrix/decomposition/qr.rs
+++ b/src/matrix/decomposition/qr.rs
@@ -226,7 +226,7 @@ impl<T> HouseholderQr<T> where T: Float {
     /// [HouseholderComposition](struct.HouseholderComposition.html)
     /// operator.
     pub fn q(&self) -> HouseholderComposition<T> {
-        householder::create_composition(&self.qr, &self.tau)
+        householder::create_composition(&self.qr, &self.tau, 0)
     }
 
     /// Computes the *thin* (or reduced) QR decomposition.
@@ -264,7 +264,7 @@ impl<T> HouseholderQr<T> where T: Float {
                 r1: qr.r
             }
         } else {
-            let composition = householder::create_composition(&self.qr, &self.tau);
+            let composition = householder::create_composition(&self.qr, &self.tau, 0);
             let q1 = composition.first_k_columns(n);
             let r1 = extract_r1(&self.qr);
             ThinQR {
@@ -291,7 +291,7 @@ impl<T: Float> Decomposition for HouseholderQr<T> {
 
 fn assemble_q<T: Float>(qr: &Matrix<T>, tau: &Vec<T>) -> Matrix<T> {
     let m = qr.rows();
-    let q_operator = householder::create_composition(qr, tau);
+    let q_operator = householder::create_composition(qr, tau, 0);
     q_operator.first_k_columns(m)
 }
 

From b02f8a073ac7f8a7431b5dd932d0ff954591c3b1 Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Fri, 5 May 2017 17:44:50 +0200
Subject: [PATCH 06/10] Implement HessenbergDecomposition

---
 src/matrix/decomposition/hessenberg.rs | 231 ++++++++++++++++++++++++-
 src/matrix/decomposition/mod.rs        |   1 +
 2 files changed, 230 insertions(+), 2 deletions(-)

diff --git a/src/matrix/decomposition/hessenberg.rs b/src/matrix/decomposition/hessenberg.rs
index b18ff35..80a8677 100644
--- a/src/matrix/decomposition/hessenberg.rs
+++ b/src/matrix/decomposition/hessenberg.rs
@@ -1,9 +1,171 @@
 use matrix::{Matrix, BaseMatrix, BaseMatrixMut, MatrixSlice, MatrixSliceMut};
+use matrix::decomposition::{Decomposition, HouseholderReflection,
+                            HouseholderComposition};
+use matrix::decomposition::householder;
 use error::{Error, ErrorKind};
+use vector::Vector;
 
 use std::any::Any;
 
-use libnum::{Float};
+use libnum::{Zero, Float};
+
+/// Result of unpacking a
+/// [Hessenberg decomposition](struct.HessenbergDecomposition.html).
+#[derive(Debug, Clone)]
+pub struct QH<T> {
+    /// The orthogonal factor `Q` in the decomposition.
+    pub q: Matrix<T>,
+    /// The upper Hessenberg factor `H` in the decomposition.
+    pub h: Matrix<T>
+}
+
+/// Upper Hessenberg decomposition of real square matrices.
+///
+/// Given any real square `m x m` matrix `A`, there exists an orthogonal
+/// `m x m` matrix `Q` and an `m x m` *upper* Hessenberg matrix `H` such that
+///
+/// ```text
+/// A = Q H Qᵀ.                                         (1)
+/// ```
+///
+/// An Upper Hessenberg matrix `H` is characterized by the property that
+/// `H[i, j] = 0` for any `i > j + 1`. For example, for a 4x4 matrix `H`, this
+/// means that `H` has a the following pattern of zeros and (possibly)
+/// non-zeros `x`:
+///
+/// ```text
+/// [ x x x x ]
+/// [ x x x x ]
+/// [ 0 x x x ]
+/// [ 0 0 x x ]
+/// ```
+///
+/// Hessenberg matrices are of particular interest because their eigenvalues
+/// and eigenvectors can be computed much more efficiently than for a general
+/// non-symmetric matrix, and because the Hessenberg matrix `H` corresponding
+/// to any matrix `A` is unitarily similar to `A`, the eigenvalues of `H`
+/// coincide with the eigenvalues of `A`.
+///
+/// As in the case of [HouseholderQr](struct.HouseholderQr.html),
+/// the `Q` factor is internally stored implicitly in terms of a sequence of
+/// Householder transformations. The implicit `Q` matrix can be accessed by calling
+/// the `.q()` method on the `HessenbergDecomposition` instance.
+///
+/// # Examples
+///
+/// ```rust
+/// # #[macro_use] extern crate rulinalg; fn main() {
+/// use rulinalg::matrix::decomposition::{Decomposition, HessenbergDecomposition, QH};
+///
+/// let x = matrix![5.0, 2.0, 0.0, -2.0;
+///                 3.0, 3.0, 4.0,  3.0;
+///                -3.0, 1.0, 2.0,  1.0;
+///                 2.0, 3.0, 2.0, -3.0];
+///
+/// // We only clone `x` here so that we can compare for equality later
+/// let QH { q, h } = HessenbergDecomposition::decompose(x.clone()).unpack();
+///
+/// assert_matrix_eq!(x, &q * h * q.transpose(), comp = float, ulp = 100);
+/// # }
+/// ```
+///
+/// # Internal storage format
+///
+/// The internal storage format is very similar to that of
+/// [HouseholderQr](struct.HouseholderQr.html). Each Householder reflector
+/// is stored compactly in the elements of the matrix which correspond
+/// to zeros in the Heissenberg factor `H`. In addition, a vector of length
+/// `n - 1` is allocated to hold the multipliers for the Householder vectors.
+#[derive(Debug, Clone)]
+pub struct HessenbergDecomposition<T> {
+    qh: Matrix<T>,
+    tau: Vec<T>
+}
+
+impl<T> HessenbergDecomposition<T> where T: Float {
+    /// Computes the upper Hessenberg decomposition
+    /// for a given square matrix.
+    ///
+    /// # Panics
+    ///
+    /// - The matrix is not square.
+    pub fn decompose(matrix: Matrix<T>) -> HessenbergDecomposition<T> {
+        assert!(matrix.rows() == matrix.cols(), "Matrix must be square for Hessenberg decomposition.");
+        let n = matrix.rows();
+
+        // n - 1 is the number of elements along the sub-1 diagonal
+        // (which characterizes the Hessenberg matrix)
+        let n1 = n.saturating_sub(1);
+
+        let mut qh = matrix;
+        let mut tau = vec![T::zero(); n1];
+
+        let mut buffer = vec![T::zero(); n1];
+        let mut multiply_buffer = vec![T::zero(); n];
+
+        for j in 0 .. n1 {
+            buffer.truncate(n1 - j);
+
+            // Compute the Householder matrix Q_j and left-apply it to the
+            // bottom right corner (we don't need to consider the columns
+            // before j, since they are implicitly zero)
+            let house = {
+                let mut bottom_right = qh.sub_slice_mut([j + 1, j], n1 - j, n - j);
+                bottom_right.col(0).clone_into_slice(&mut buffer);
+
+                let house = HouseholderReflection::compute(Vector::new(buffer));
+                house.buffered_left_multiply_into(&mut bottom_right, &mut multiply_buffer);
+                house.store_in_col(&mut bottom_right.col_mut(0));
+                house
+            };
+
+            // Apply Q_j to the trailing columns
+            let mut trailing_cols = qh.sub_slice_mut([0, j + 1], n, n1 - j);
+            house.buffered_right_multiply_into(&mut trailing_cols, &mut multiply_buffer);
+
+            tau[j] = house.tau();
+            buffer = house.into_vector().into_vec();
+
+        }
+
+        HessenbergDecomposition {
+            qh: qh,
+            tau: tau
+        }
+    }
+
+    /// Returns the orthogonal factor Q as an instance
+    /// of a [HouseholderComposition](struct.HouseholderComposition.html)
+    /// operator.
+    pub fn q(&self) -> HouseholderComposition<T> {
+        householder::create_composition(&self.qh, &self.tau, 1)
+    }
+}
+
+impl<T> Decomposition for HessenbergDecomposition<T> where T: Float {
+    type Factors = QH<T>;
+
+    fn unpack(self) -> Self::Factors {
+        let q = self.q().first_k_columns(self.qh.rows());
+        let mut h = self.qh;
+        nullify_lower_portion(&mut h);
+
+        QH {
+            q: q,
+            h: h
+        }
+    }
+}
+
+/// Makes the portion of the matrix which is not part of the
+/// Upper Hessenberg part identically zero.
+fn nullify_lower_portion<T>(matrix: &mut Matrix<T>) where T: Zero {
+    for (i, mut row) in matrix.row_iter_mut().enumerate().skip(1) {
+        for element in row.raw_slice_mut().iter_mut().take(i - 1) {
+            *element = T::zero();
+        }
+    }
+}
 
 impl<T: Any + Float> Matrix<T> {
     /// Returns H, where H is the upper hessenberg form.
@@ -144,7 +306,12 @@ impl<T: Any + Float> Matrix<T> {
 
 #[cfg(test)]
 mod tests {
-    use matrix::Matrix;
+    use super::QH;
+    use super::HessenbergDecomposition;
+
+    use matrix::{BaseMatrix, Matrix};
+    use matrix::decomposition::Decomposition;
+    use testsupport::is_upper_hessenberg;
 
     #[test]
     #[should_panic]
@@ -161,4 +328,64 @@ mod tests {
 
         let _ = a.upper_hess_decomp();
     }
+
+    fn verify_hessenberg(x: Matrix<f64>) {
+        let QH {ref h, ref q } = HessenbergDecomposition::decompose(x.clone()).unpack();
+
+        assert!(is_upper_hessenberg(h));
+
+        let identity = Matrix::identity(h.rows());
+
+        // Orthogonality. For simplicity, we use a fixed absolute tolerance
+        // and expect elements to be in the vicinity of 1.0
+        assert_matrix_eq!(q.transpose() * q, identity, comp = abs, tol = 1e-14);
+        assert_matrix_eq!(q * q.transpose(), identity, comp = abs, tol = 1e-14);
+
+        // Reconstruction
+        assert_matrix_eq!(q * h * q.transpose(), x, comp = abs, tol = 1e-14);
+    }
+
+    #[test]
+    fn hessenberg_decomposition() {
+        {
+            let x: Matrix<f64> = matrix![];
+            verify_hessenberg(x);
+        }
+
+        {
+            let x = matrix![3.0];
+            verify_hessenberg(x);
+        }
+
+        {
+            let x = matrix![3.0,  2.0;
+                            5.0,  1.0];
+            verify_hessenberg(x);
+        }
+
+        {
+            let x = matrix![3.0,  2.0, -5.0;
+                            5.0,  1.0,  2.0;
+                            3.0, -2.0,  0.0];
+            verify_hessenberg(x);
+        }
+
+        {
+            let x = matrix![3.0,  2.0, -5.0,  2.0;
+                            5.0,  1.0,  2.0, -2.0;
+                            3.0, -2.0,  0.0,  0.0;
+                            4.0,  3.0, -1.0,  3.0];
+            verify_hessenberg(x);
+        }
+
+        {
+            let x = matrix![1.0,  3.0, -5.0,  2.0, 2.0;
+                            0.0,  0.0,  5.0, -1.0, 4.0;
+                            3.0, -4.0,  0.0,  0.0, 5.0;
+                            0.0,  3.0, -1.0,  3.0, 1.0;
+                            2.0, -4.0,  3.0,  2.0, 3.0];
+            verify_hessenberg(x);
+        }
+    }
+
 }
diff --git a/src/matrix/decomposition/mod.rs b/src/matrix/decomposition/mod.rs
index d2f9d93..ecfc59e 100644
--- a/src/matrix/decomposition/mod.rs
+++ b/src/matrix/decomposition/mod.rs
@@ -147,6 +147,7 @@ pub use self::householder::HouseholderComposition;
 pub use self::lu::{PartialPivLu, LUP, FullPivLu, LUPQ};
 pub use self::cholesky::Cholesky;
 pub use self::qr::{HouseholderQr, QR, ThinQR};
+pub use self::hessenberg::{HessenbergDecomposition, QH};
 
 use libnum::{Float};
 

From 029755e44f3e906bed969609d9efa6e82fe1d5ea Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Fri, 5 May 2017 18:23:26 +0200
Subject: [PATCH 07/10] Benchmarks for Hessenberg decomposition

---
 benches/lib.rs               |  1 +
 benches/linalg/hessenberg.rs | 73 ++++++++++++++++++++++++++++++++++++
 2 files changed, 74 insertions(+)
 create mode 100644 benches/linalg/hessenberg.rs

diff --git a/benches/lib.rs b/benches/lib.rs
index f693a75..d2a6e11 100644
--- a/benches/lib.rs
+++ b/benches/lib.rs
@@ -16,5 +16,6 @@ pub mod linalg {
 	mod norm;
 	mod triangular;
 	mod permutation;
+	mod hessenberg;
 	pub mod util;
 }
diff --git a/benches/linalg/hessenberg.rs b/benches/linalg/hessenberg.rs
new file mode 100644
index 0000000..3086f91
--- /dev/null
+++ b/benches/linalg/hessenberg.rs
@@ -0,0 +1,73 @@
+use rulinalg::matrix::decomposition::{Decomposition, HessenbergDecomposition};
+
+use linalg::util::reproducible_random_matrix;
+
+use test::Bencher;
+
+#[bench]
+fn hessenberg_decomposition_decompose_100x100(b: &mut Bencher) {
+    let x = reproducible_random_matrix(100, 100);
+    b.iter(|| {
+        HessenbergDecomposition::decompose(x.clone())
+    })
+}
+
+#[bench]
+fn hessenberg_decomposition_decompose_unpack_100x100(b: &mut Bencher) {
+    let x = reproducible_random_matrix(100, 100);
+    b.iter(|| {
+        HessenbergDecomposition::decompose(x.clone()).unpack()
+    })
+}
+
+#[bench]
+fn hessenberg_decomposition_decompose_200x200(b: &mut Bencher) {
+    let x = reproducible_random_matrix(200, 200);
+    b.iter(|| {
+        HessenbergDecomposition::decompose(x.clone())
+    })
+}
+
+#[bench]
+fn hessenberg_decomposition_decompose_unpack_200x200(b: &mut Bencher) {
+    let x = reproducible_random_matrix(200, 200);
+    b.iter(|| {
+        HessenbergDecomposition::decompose(x.clone()).unpack()
+    })
+}
+
+#[bench]
+#[allow(deprecated)]
+fn upper_hessenberg_100x100(b: &mut Bencher) {
+    let x = reproducible_random_matrix(100, 100);
+    b.iter(|| {
+        x.clone().upper_hessenberg()
+    })
+}
+
+#[bench]
+#[allow(deprecated)]
+fn upper_hess_decomp_100x100(b: &mut Bencher) {
+    let x = reproducible_random_matrix(100, 100);
+    b.iter(|| {
+        x.clone().upper_hess_decomp()
+    })
+}
+
+#[bench]
+#[allow(deprecated)]
+fn upper_hessenberg_200x200(b: &mut Bencher) {
+    let x = reproducible_random_matrix(200, 200);
+    b.iter(|| {
+        x.clone().upper_hessenberg()
+    })
+}
+
+#[bench]
+#[allow(deprecated)]
+fn upper_hess_decomp_200x200(b: &mut Bencher) {
+    let x = reproducible_random_matrix(200, 200);
+    b.iter(|| {
+        x.clone().upper_hess_decomp()
+    })
+}

From afe50acaf139fdbb28a2617551d5db46543ac3bf Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Sun, 7 May 2017 11:23:01 +0200
Subject: [PATCH 08/10] Mention Hessenberg in docs for HouseholderComposition

---
 src/matrix/decomposition/householder.rs | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/matrix/decomposition/householder.rs b/src/matrix/decomposition/householder.rs
index 0458555..735454c 100644
--- a/src/matrix/decomposition/householder.rs
+++ b/src/matrix/decomposition/householder.rs
@@ -201,8 +201,10 @@ impl<T: Float> HouseholderReflection<T> {
 /// Q = Q_1 * Q_2 * ... * Q_p
 /// ```
 ///
-/// as explained in the documentation for
-/// [HouseholderQr](struct.HouseholderQr.html).
+/// where `p` is the number of transformations. See the
+/// documentation for [HouseholderQr](struct.HouseholderQr.html)
+/// and [HessenbergDecomposition](struct.HessenbergDecomposition.html)
+/// for more details.
 #[derive(Debug, Clone)]
 pub struct HouseholderComposition<'a, T> where T: 'a {
     storage: &'a Matrix<T>,

From aba73af236254af77e2fcdba22d3009a30c18cc8 Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Sun, 7 May 2017 11:37:08 +0200
Subject: [PATCH 09/10] Add HessenbergDecomposition to table of decompositions

---
 src/matrix/decomposition/mod.rs | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/src/matrix/decomposition/mod.rs b/src/matrix/decomposition/mod.rs
index ecfc59e..e737c4c 100644
--- a/src/matrix/decomposition/mod.rs
+++ b/src/matrix/decomposition/mod.rs
@@ -110,6 +110,12 @@
 //! <td></td>
 //! </tr>
 //!
+//! <tr>
+//! <td><a href="struct.HessenbergDecomposition.html">HessenbergDecomposition</a></td>
+//! <td>Square</td>
+//! <td></td>
+//! </tr>
+//!
 //! </tbody>
 //! </table>
 

From bcbd73fdf6f7cba4b5b6a763ed283e5f029b5820 Mon Sep 17 00:00:00 2001
From: Andreas Longva <andreas.b.longva@gmail.com>
Date: Sun, 7 May 2017 11:43:23 +0200
Subject: [PATCH 10/10] Deprecate old Hessenberg functions

---
 src/matrix/decomposition/eigen.rs      |  2 ++
 src/matrix/decomposition/hessenberg.rs | 15 +++++++++++++++
 2 files changed, 17 insertions(+)

diff --git a/src/matrix/decomposition/eigen.rs b/src/matrix/decomposition/eigen.rs
index eb328f2..b6466e6 100644
--- a/src/matrix/decomposition/eigen.rs
+++ b/src/matrix/decomposition/eigen.rs
@@ -82,6 +82,7 @@ impl<T: Any + Float + Signed> Matrix<T> {
                       "Francis shift only works on matrices greater than 2x2.");
         debug_assert!(n == self.cols, "Matrix must be square for Francis shift.");
 
+        #[allow(deprecated)]
         let mut h = try!(self
             .upper_hessenberg()
             .map_err(|_| Error::new(ErrorKind::DecompFailure, "Could not compute eigenvalues.")));
@@ -223,6 +224,7 @@ impl<T: Any + Float + Signed> Matrix<T> {
                       "Francis shift only works on matrices greater than 2x2.");
         debug_assert!(n == self.cols, "Matrix must be square for Francis shift.");
 
+        #[allow(deprecated)]
         let (u, mut h) = try!(self.upper_hess_decomp().map_err(|_| {
             Error::new(ErrorKind::DecompFailure,
                        "Could not compute eigen decomposition.")
diff --git a/src/matrix/decomposition/hessenberg.rs b/src/matrix/decomposition/hessenberg.rs
index 80a8677..ef1366e 100644
--- a/src/matrix/decomposition/hessenberg.rs
+++ b/src/matrix/decomposition/hessenberg.rs
@@ -173,6 +173,11 @@ impl<T: Any + Float> Matrix<T> {
     /// If the transformation matrix is also required, you should
     /// use `upper_hess_decomp`.
     ///
+    /// Note: This function is deprecated and will be removed in a future
+    /// release. See
+    /// [HessenbergDecomposition](decomposition/struct.HessenbergDecomposition.html)
+    /// for its replacement.
+    ///
     /// # Examples
     ///
     /// ```
@@ -196,6 +201,7 @@ impl<T: Any + Float> Matrix<T> {
     /// # Failures
     ///
     /// - The matrix cannot be reduced to upper hessenberg form.
+    #[deprecated]
     pub fn upper_hessenberg(mut self) -> Result<Matrix<T>, Error> {
         let n = self.rows;
         assert!(n == self.cols,
@@ -249,6 +255,11 @@ impl<T: Any + Float> Matrix<T> {
     ///
     /// Note: The current transform matrix seems broken...
     ///
+    /// Note: This function is deprecated and will be removed in a future
+    /// release. See
+    /// [HessenbergDecomposition](decomposition/struct.HessenbergDecomposition.html)
+    /// for its replacement.
+    ///
     /// # Examples
     ///
     /// ```
@@ -273,6 +284,7 @@ impl<T: Any + Float> Matrix<T> {
     /// # Failures
     ///
     /// - The matrix cannot be reduced to upper hessenberg form.
+    #[deprecated]
     pub fn upper_hess_decomp(self) -> Result<(Matrix<T>, Matrix<T>), Error> {
         let n = self.rows;
         assert!(n == self.cols,
@@ -300,6 +312,7 @@ impl<T: Any + Float> Matrix<T> {
         }
 
         // Now we reduce to upper hessenberg
+        #[allow(deprecated)]
         Ok((transform, try!(self.upper_hessenberg())))
     }
 }
@@ -315,6 +328,7 @@ mod tests {
 
     #[test]
     #[should_panic]
+    #[allow(deprecated)]
     fn test_non_square_upper_hessenberg() {
         let a: Matrix<f64> = Matrix::ones(2, 3);
 
@@ -323,6 +337,7 @@ mod tests {
 
     #[test]
     #[should_panic]
+    #[allow(deprecated)]
     fn test_non_square_upper_hess_decomp() {
         let a: Matrix<f64> = Matrix::ones(2, 3);