diff --git a/src/householder.rs b/src/householder.rs
new file mode 100644
index 00000000..6ef0bd2e
--- /dev/null
+++ b/src/householder.rs
@@ -0,0 +1,84 @@
+use crate::{inner::*, norm::*, types::*};
+use ndarray::*;
+
+/// Iterative orthogonalizer using Householder reflection
+#[derive(Debug, Clone)]
+pub struct Householder<A> {
+    dim: usize,
+    v: Vec<Array1<A>>,
+}
+
+impl<A: Scalar> Householder<A> {
+    pub fn new(dim: usize) -> Self {
+        Householder { dim, v: Vec::new() }
+    }
+
+    pub fn len(&self) -> usize {
+        self.v.len()
+    }
+
+    /// Take a Reflection `P = I - 2ww^T`
+    fn reflect<S: DataMut<Elem = A>>(&self, k: usize, a: &mut ArrayBase<S, Ix1>) {
+        assert!(k < self.v.len());
+        assert_eq!(a.len(), self.dim);
+        let w = self.v[k].slice(s![k..]);
+        let c = A::from(2.0).unwrap() * w.inner(&a.slice(s![k..]));
+        for l in k..self.dim {
+            a[l] -= c * w[l];
+        }
+    }
+
+    /// Orghotonalize input vector by reflectors
+    ///
+    /// Panic
+    /// -------
+    /// - if the size of the input array mismaches to the dimension
+    pub fn orthogonalize<S>(&self, a: &mut ArrayBase<S, Ix1>) -> A::Real
+    where
+        A: Lapack,
+        S: DataMut<Elem = A>,
+    {
+        for k in 0..self.len() {
+            self.reflect(k, a);
+        }
+        // residual norm
+        a.slice(s![self.len()..]).norm_l2()
+    }
+
+    /// Orghotonalize input vector by reflectors
+    ///
+    /// ```rust
+    /// # use ndarray::*;
+    /// # use ndarray_linalg::*;
+    /// let mut mgs = arnoldi::MGS::new(3);
+    /// let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap();
+    /// close_l2(&coef, &array![1.0], 1e-9).unwrap();
+    ///
+    /// let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap();
+    /// close_l2(&coef, &array![1.0, 1.0], 1e-9).unwrap();
+    ///
+    /// assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err());  // Fail if the vector is linearly dependend
+    ///
+    /// if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) {
+    ///     close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9).unwrap(); // You can get coefficients of dependent vector
+    /// }
+    /// ```
+    ///
+    /// Panic
+    /// -------
+    /// - if the size of the input array mismaches to the dimension
+    ///
+    pub fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    where
+        A: Lapack,
+        S: DataMut<Elem = A>,
+    {
+        let residual = self.orthogonalize(&mut a);
+        let coef = a.slice(s![..self.len()]).into_owned();
+        if residual < rtol {
+            return Err(coef);
+        }
+        self.v.push(a.into_owned());
+        Ok(coef)
+    }
+}
diff --git a/src/inner.rs b/src/inner.rs
new file mode 100644
index 00000000..299e9669
--- /dev/null
+++ b/src/inner.rs
@@ -0,0 +1,21 @@
+use crate::types::*;
+use ndarray::*;
+
+pub trait Inner<S: Data> {
+    /// Inner product `(self.conjugate, rhs)`
+    fn inner(&self, rhs: &ArrayBase<S, Ix1>) -> S::Elem;
+}
+
+impl<A, S1, S2> Inner<S1> for ArrayBase<S2, Ix1>
+where
+    A: Scalar,
+    S1: Data<Elem = A>,
+    S2: Data<Elem = A>,
+{
+    fn inner(&self, rhs: &ArrayBase<S1, Ix1>) -> A {
+        Zip::from(self)
+            .and(rhs)
+            .fold_while(A::zero(), |acc, s, r| FoldWhile::Continue(acc + s.conj() * *r))
+            .into_inner()
+    }
+}
diff --git a/src/krylov/householder.rs b/src/krylov/householder.rs
new file mode 100644
index 00000000..9e7cef9d
--- /dev/null
+++ b/src/krylov/householder.rs
@@ -0,0 +1,107 @@
+use super::*;
+use crate::{inner::*, norm::*};
+use num_traits::Zero;
+
+/// Iterative orthogonalizer using Householder reflection
+#[derive(Debug, Clone)]
+pub struct Householder<A: Scalar> {
+    /// Dimension of orthogonalizer
+    dim: usize,
+
+    /// Store Householder reflector.
+    ///
+    /// The coefficient is copied into another array, and this does not contain
+    v: Vec<Array1<A>>,
+}
+
+impl<A: Scalar> Householder<A> {
+    /// Create a new orthogonalizer
+    pub fn new(dim: usize) -> Self {
+        Householder { dim, v: Vec::new() }
+    }
+
+    /// Take a Reflection `P = I - 2ww^T`
+    fn reflect<S: DataMut<Elem = A>>(&self, k: usize, a: &mut ArrayBase<S, Ix1>) {
+        assert!(k < self.v.len());
+        assert_eq!(a.len(), self.dim, "Input array size mismaches to the dimension");
+
+        let w = self.v[k].slice(s![k..]);
+        let mut a_slice = a.slice_mut(s![k..]);
+        let c = A::from(2.0).unwrap() * w.inner(&a_slice);
+        for l in 0..self.dim - k {
+            a_slice[l] -= c * w[l];
+        }
+    }
+}
+
+impl<A: Scalar + Lapack> Orthogonalizer for Householder<A> {
+    type Elem = A;
+
+    fn dim(&self) -> usize {
+        self.dim
+    }
+
+    fn len(&self) -> usize {
+        self.v.len()
+    }
+
+    fn orthogonalize<S>(&self, a: &mut ArrayBase<S, Ix1>) -> A::Real
+    where
+        S: DataMut<Elem = A>,
+    {
+        for k in 0..self.len() {
+            self.reflect(k, a);
+        }
+        if self.is_full() {
+            return Zero::zero();
+        }
+        // residual norm
+        a.slice(s![self.len()..]).norm_l2()
+    }
+
+    fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    where
+        S: DataMut<Elem = A>,
+    {
+        assert_eq!(a.len(), self.dim);
+        let k = self.len();
+        let alpha = self.orthogonalize(&mut a);
+        let mut coef = Array::zeros(k + 1);
+        for i in 0..k {
+            coef[i] = a[i];
+        }
+        if alpha < rtol {
+            // linearly dependent
+            coef[k] = A::from_real(alpha);
+            return Err(coef);
+        }
+
+        // Add reflector
+        assert!(k < a.len()); // this must hold because `alpha == 0` if k >= a.len()
+        let alpha = if a[k].abs() > Zero::zero() {
+            a[k].mul_real(alpha / a[k].abs())
+        } else {
+            A::from_real(alpha)
+        };
+        coef[k] = alpha;
+
+        a[k] -= alpha;
+        let norm = a.slice(s![k..]).norm_l2();
+        azip!(mut a (a.slice_mut(s![..k])) in { *a = Zero::zero() }); // this can be omitted
+        azip!(mut a (a.slice_mut(s![k..])) in { *a = a.div_real(norm) });
+        self.v.push(a.into_owned());
+        Ok(coef)
+    }
+
+    fn get_q(&self) -> Q<A> {
+        assert!(self.len() > 0);
+        let mut a = Array::zeros((self.dim(), self.len()));
+        for (i, mut col) in a.axis_iter_mut(Axis(1)).enumerate() {
+            col[i] = A::one();
+            for l in (0..self.len()).rev() {
+                self.reflect(l, &mut col);
+            }
+        }
+        a
+    }
+}
diff --git a/src/krylov/mgs.rs b/src/krylov/mgs.rs
new file mode 100644
index 00000000..ae9aa1a6
--- /dev/null
+++ b/src/krylov/mgs.rs
@@ -0,0 +1,101 @@
+use super::*;
+use crate::{generate::*, inner::*, norm::Norm};
+
+/// Iterative orthogonalizer using modified Gram-Schmit procedure
+#[derive(Debug, Clone)]
+pub struct MGS<A> {
+    /// Dimension of base space
+    dim: usize,
+    /// Basis of spanned space
+    q: Vec<Array1<A>>,
+}
+
+impl<A: Scalar> MGS<A> {
+    pub fn new(dim: usize) -> Self {
+        Self { dim, q: Vec::new() }
+    }
+
+    fn ortho<S>(&self, a: &mut ArrayBase<S, Ix1>) -> Array1<A>
+    where
+        A: Lapack,
+        S: DataMut<Elem = A>,
+    {
+        assert_eq!(a.len(), self.dim);
+        let mut coef = Array1::zeros(self.q.len() + 1);
+        for i in 0..self.q.len() {
+            let q = &self.q[i];
+            let c = q.inner(&a);
+            azip!(mut a (&mut *a), q (q) in { *a = *a - c * q } );
+            coef[i] = c;
+        }
+        let nrm = a.norm_l2();
+        coef[self.q.len()] = A::from(nrm).unwrap();
+        coef
+    }
+}
+
+impl<A: Scalar + Lapack> Orthogonalizer for MGS<A> {
+    type Elem = A;
+
+    fn dim(&self) -> usize {
+        self.dim
+    }
+
+    fn len(&self) -> usize {
+        self.q.len()
+    }
+
+    fn orthogonalize<S>(&self, a: &mut ArrayBase<S, Ix1>) -> A::Real
+    where
+        S: DataMut<Elem = A>,
+    {
+        let coef = self.ortho(a);
+        // Write coefficients into `a`
+        azip!(mut a (a.slice_mut(s![0..self.len()])), coef in { *a = coef });
+        // 0-fill for remaining
+        azip!(mut a (a.slice_mut(s![self.len()..])) in { *a = A::zero() });
+        coef[self.len()].re()
+    }
+
+    fn append<S>(&mut self, mut a: ArrayBase<S, Ix1>, rtol: A::Real) -> Result<Array1<A>, Array1<A>>
+    where
+        S: DataMut<Elem = A>,
+    {
+        let coef = self.ortho(&mut a);
+        let nrm = coef[self.len()].re();
+        if nrm < rtol {
+            // Linearly dependent
+            return Err(coef);
+        }
+        azip!(mut a in { *a = *a / A::from_real(nrm) });
+        self.q.push(a.into_owned());
+        Ok(coef)
+    }
+
+    fn get_q(&self) -> Q<A> {
+        hstack(&self.q).unwrap()
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+    use crate::assert::*;
+
+    #[test]
+    fn mgs_append() {
+        let mut mgs = MGS::new(3);
+        let coef = mgs.append(array![0.0, 1.0, 0.0], 1e-9).unwrap();
+        close_l2(&coef, &array![1.0], 1e-9);
+
+        let coef = mgs.append(array![1.0, 1.0, 0.0], 1e-9).unwrap();
+        close_l2(&coef, &array![1.0, 1.0], 1e-9);
+
+        assert!(mgs.append(array![1.0, 2.0, 0.0], 1e-9).is_err());
+
+        if let Err(coef) = mgs.append(array![1.0, 2.0, 0.0], 1e-9) {
+            close_l2(&coef, &array![2.0, 1.0, 0.0], 1e-9);
+        }
+    }
+
+}
diff --git a/src/krylov/mod.rs b/src/krylov/mod.rs
new file mode 100644
index 00000000..fb5b15a4
--- /dev/null
+++ b/src/krylov/mod.rs
@@ -0,0 +1,156 @@
+use crate::types::*;
+use ndarray::*;
+
+mod householder;
+mod mgs;
+
+pub use householder::Householder;
+pub use mgs::MGS;
+
+/// Q-matrix (unitary)
+pub type Q<A> = Array2<A>;
+/// R-matrix (upper triangle)
+pub type R<A> = Array2<A>;
+
+pub trait Orthogonalizer {
+    type Elem: Scalar;
+
+    fn dim(&self) -> usize;
+    fn len(&self) -> usize;
+
+    /// check if the basis spans entire space
+    fn is_full(&self) -> bool {
+        self.len() == self.dim()
+    }
+
+    fn is_empty(&self) -> bool {
+        self.len() == 0
+    }
+
+    /// Orthogonalize given vector
+    ///
+    /// Returns
+    /// --------
+    /// Residual norm
+    ///
+    /// Panic
+    /// -------
+    /// - if the size of the input array mismaches to the dimension
+    ///
+    fn orthogonalize<S>(&self, a: &mut ArrayBase<S, Ix1>) -> <Self::Elem as Scalar>::Real
+    where
+        S: DataMut<Elem = Self::Elem>;
+
+    /// Add new vector if the residual is larger than relative tolerance
+    ///
+    /// Returns
+    /// --------
+    /// Coefficients to the `i`-th Q-vector
+    ///
+    /// - The size of array must be `self.len() + 1`
+    /// - The last element is the residual norm of input vector
+    ///
+    /// Panic
+    /// -------
+    /// - if the size of the input array mismaches to the dimension
+    ///
+    fn append<S>(
+        &mut self,
+        a: ArrayBase<S, Ix1>,
+        rtol: <Self::Elem as Scalar>::Real,
+    ) -> Result<Array1<Self::Elem>, Array1<Self::Elem>>
+    where
+        S: DataMut<Elem = Self::Elem>;
+
+    /// Get Q-matrix of generated basis
+    fn get_q(&self) -> Q<Self::Elem>;
+}
+
+/// Strategy for linearly dependent vectors appearing in iterative QR decomposition
+#[derive(Clone, Copy, Debug, Eq, PartialEq)]
+pub enum Strategy {
+    /// Terminate iteration if dependent vector comes
+    Terminate,
+
+    /// Skip dependent vector
+    Skip,
+
+    /// Orghotonalize dependent vector without adding to Q,
+    /// thus R must be non-regular like following:
+    ///
+    /// ```ignore
+    /// x x x x x
+    /// 0 x x x x
+    /// 0 0 0 x x
+    /// 0 0 0 0 x
+    /// 0 0 0 0 0   // 0-filled to be square matrix
+    /// ```
+    Full,
+}
+
+/// Online QR decomposition using arbitary orthogonalizer
+pub fn qr<A, S>(
+    iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
+    mut ortho: impl Orthogonalizer<Elem = A>,
+    rtol: A::Real,
+    strategy: Strategy,
+) -> (Q<A>, R<A>)
+where
+    A: Scalar + Lapack,
+    S: Data<Elem = A>,
+{
+    assert_eq!(ortho.len(), 0);
+
+    let mut coefs = Vec::new();
+    for a in iter {
+        match ortho.append(a.into_owned(), rtol) {
+            Ok(coef) => coefs.push(coef),
+            Err(coef) => match strategy {
+                Strategy::Terminate => break,
+                Strategy::Skip => continue,
+                Strategy::Full => coefs.push(coef),
+            },
+        }
+    }
+    let n = ortho.len();
+    let m = coefs.len();
+    let mut r = Array2::zeros((n, m).f());
+    for j in 0..m {
+        for i in 0..n {
+            if i < coefs[j].len() {
+                r[(i, j)] = coefs[j][i];
+            }
+        }
+    }
+    (ortho.get_q(), r)
+}
+
+/// Online QR decomposition using modified Gram-Schmit
+pub fn mgs<A, S>(
+    iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
+    dim: usize,
+    rtol: A::Real,
+    strategy: Strategy,
+) -> (Q<A>, R<A>)
+where
+    A: Scalar + Lapack,
+    S: Data<Elem = A>,
+{
+    let mgs = MGS::new(dim);
+    qr(iter, mgs, rtol, strategy)
+}
+
+/// Online QR decomposition using modified Gram-Schmit
+pub fn householder<A, S>(
+    iter: impl Iterator<Item = ArrayBase<S, Ix1>>,
+    dim: usize,
+    rtol: A::Real,
+    strategy: Strategy,
+) -> (Q<A>, R<A>)
+where
+    A: Scalar + Lapack,
+    S: Data<Elem = A>,
+{
+    let h = Householder::new(dim);
+    qr(iter, h, rtol, strategy)
+}
diff --git a/src/lib.rs b/src/lib.rs
index 1a55a082..ed4d910d 100644
--- a/src/lib.rs
+++ b/src/lib.rs
@@ -43,6 +43,8 @@ pub mod diagonal;
 pub mod eigh;
 pub mod error;
 pub mod generate;
+pub mod inner;
+pub mod krylov;
 pub mod lapack;
 pub mod layout;
 pub mod norm;
diff --git a/tests/householder.rs b/tests/householder.rs
new file mode 100644
index 00000000..adc4d9b2
--- /dev/null
+++ b/tests/householder.rs
@@ -0,0 +1,96 @@
+use ndarray::*;
+use ndarray_linalg::{krylov::*, *};
+
+fn over<A: Scalar + Lapack>(rtol: A::Real) {
+    const N: usize = 4;
+    let a: Array2<A> = random((N, N * 2));
+
+    // Terminate
+    let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
+    let a_sub = a.slice(s![.., 0..N]);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I");
+    assert_close_l2!(&q.dot(&r), &a_sub, rtol; "Check A = QR");
+
+    // Skip
+    let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip);
+    let a_sub = a.slice(s![.., 0..N]);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
+    assert_close_l2!(&q.dot(&r), &a_sub, rtol);
+
+    // Full
+    let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Full);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
+    assert_close_l2!(&q.dot(&r), &a, rtol);
+}
+
+#[test]
+fn over_f32() {
+    over::<f32>(1e-5);
+}
+#[test]
+fn over_f64() {
+    over::<f64>(1e-9);
+}
+#[test]
+fn over_c32() {
+    over::<c32>(1e-5);
+}
+#[test]
+fn over_c64() {
+    over::<c64>(1e-9);
+}
+
+fn full<A: Scalar + Lapack>(rtol: A::Real) {
+    const N: usize = 5;
+    let a: Array2<A> = random((N, N));
+    let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol; "Check Q^H Q = I");
+    assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR");
+}
+
+#[test]
+fn full_f32() {
+    full::<f32>(1e-5);
+}
+#[test]
+fn full_f64() {
+    full::<f64>(1e-9);
+}
+#[test]
+fn full_c32() {
+    full::<c32>(1e-5);
+}
+#[test]
+fn full_c64() {
+    full::<c64>(1e-9);
+}
+
+fn half<A: Scalar + Lapack>(rtol: A::Real) {
+    const N: usize = 4;
+    let a: Array2<A> = random((N, N / 2));
+    let (q, r) = householder(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol; "Check Q^H Q = I");
+    assert_close_l2!(&q.dot(&r), &a, rtol; "Check A = QR");
+}
+
+#[test]
+fn half_f32() {
+    half::<f32>(1e-5);
+}
+#[test]
+fn half_f64() {
+    half::<f64>(1e-9);
+}
+#[test]
+fn half_c32() {
+    half::<c32>(1e-5);
+}
+#[test]
+fn half_c64() {
+    half::<c64>(1e-9);
+}
diff --git a/tests/mgs.rs b/tests/mgs.rs
new file mode 100644
index 00000000..f80bbf5e
--- /dev/null
+++ b/tests/mgs.rs
@@ -0,0 +1,96 @@
+use ndarray::*;
+use ndarray_linalg::{krylov::*, *};
+
+fn full<A: Scalar + Lapack>(rtol: A::Real) {
+    const N: usize = 5;
+    let a: Array2<A> = random((N, N));
+    let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
+    assert_close_l2!(&q.dot(&r), &a, rtol);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
+}
+
+#[test]
+fn full_f32() {
+    full::<f32>(1e-5);
+}
+#[test]
+fn full_f64() {
+    full::<f64>(1e-9);
+}
+#[test]
+fn full_c32() {
+    full::<c32>(1e-5);
+}
+#[test]
+fn full_c64() {
+    full::<c64>(1e-9);
+}
+
+fn half<A: Scalar + Lapack>(rtol: A::Real) {
+    const N: usize = 4;
+    let a: Array2<A> = random((N, N / 2));
+    let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
+    assert_close_l2!(&q.dot(&r), &a, rtol);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N / 2), rtol);
+}
+
+#[test]
+fn half_f32() {
+    half::<f32>(1e-5);
+}
+#[test]
+fn half_f64() {
+    half::<f64>(1e-9);
+}
+#[test]
+fn half_c32() {
+    half::<c32>(1e-5);
+}
+#[test]
+fn half_c64() {
+    half::<c64>(1e-9);
+}
+
+fn over<A: Scalar + Lapack>(rtol: A::Real) {
+    const N: usize = 4;
+    let a: Array2<A> = random((N, N * 2));
+
+    // Terminate
+    let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Terminate);
+    let a_sub = a.slice(s![.., 0..N]);
+    assert_close_l2!(&q.dot(&r), &a_sub, rtol);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
+
+    // Skip
+    let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Skip);
+    let a_sub = a.slice(s![.., 0..N]);
+    assert_close_l2!(&q.dot(&r), &a_sub, rtol);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
+
+    // Full
+    let (q, r) = mgs(a.axis_iter(Axis(1)), N, rtol, Strategy::Full);
+    assert_close_l2!(&q.dot(&r), &a, rtol);
+    let qc: Array2<A> = conjugate(&q);
+    assert_close_l2!(&qc.dot(&q), &Array::eye(N), rtol);
+}
+
+#[test]
+fn over_f32() {
+    over::<f32>(1e-5);
+}
+#[test]
+fn over_f64() {
+    over::<f64>(1e-9);
+}
+#[test]
+fn over_c32() {
+    over::<c32>(1e-5);
+}
+#[test]
+fn over_c64() {
+    over::<c64>(1e-9);
+}