diff --git a/Project.toml b/Project.toml
index 320b2546..2a2ffe0d 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "RobustAndOptimalControl"
 uuid = "21fd56a4-db03-40ee-82ee-a87907bee541"
 authors = ["Fredrik Bagge Carlson", "Marcus Greiff"]
-version = "0.3.4"
+version = "0.4.0"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
diff --git a/src/RobustAndOptimalControl.jl b/src/RobustAndOptimalControl.jl
index d15b9d56..e9a80064 100644
--- a/src/RobustAndOptimalControl.jl
+++ b/src/RobustAndOptimalControl.jl
@@ -24,7 +24,7 @@ using ChainRulesCore
 export show_construction, vec2sys
 include("utils.jl")
 
-export dss, hinfnorm2, h2norm, hankelnorm, nugap, νgap, baltrunc2, stab_unstab, baltrunc_unstab, baltrunc_coprime
+export dss, hinfnorm2, linfnorm2, h2norm, hankelnorm, nugap, νgap, baltrunc2, stab_unstab, baltrunc_unstab, baltrunc_coprime
 include("descriptor.jl")
 
 export ExtendedStateSpace, system_mapping, performance_mapping, noise_mapping, ssdata_e, partition, ss
diff --git a/src/canonical.jl b/src/canonical.jl
index 0242079d..a5ee3e33 100644
--- a/src/canonical.jl
+++ b/src/canonical.jl
@@ -1,4 +1,8 @@
-blockdiagonalize(A::AbstractMatrix) = cdf2rdf(eigen(A, sortby=eigsortby))
+function blockdiagonalize(A::AbstractMatrix)
+    E = eigen(A, sortby=eigsortby)
+    Db,Vb = cdf2rdf(E)
+    Db,Vb,E
+end
 
 eigsortby(λ::Real) = λ
 eigsortby(λ::Complex) = (abs(imag(λ)),real(λ))
@@ -54,7 +58,7 @@ function cdf2rdf(E::Eigen)
 end
 
 """
-    sysm, T = modal_form(sys; C1 = false)
+    sysm, T, E = modal_form(sys; C1 = false)
 
 Bring `sys` to modal form.
 
@@ -65,10 +69,12 @@ sysm ≈ similarity_transform(sys, T)
 
 If `C1`, then an additional convention for SISO systems is used, that the `C`-matrix coefficient of real eigenvalues is 1. If `C1 = false`, the `B` and `C` coefficients are chosen in a balanced fashion.
 
+`E` is an eigen factorization of `A`.
+
 See also [`hess_form`](@ref) and [`schur_form`](@ref)
 """
 function modal_form(sys; C1 = false)
-    Ab,T = blockdiagonalize(sys.A)
+    Ab,T,E = blockdiagonalize(sys.A)
     # Calling similarity_transform looks like a detour, but this implementation allows modal_form to work with any AbstractStateSpace which implements a custom method for similarity transform
     sysm = similarity_transform(sys, T)
     sysm.A .= Ab # sysm.A should already be Ab after similarity_transform, but Ab has less numerical noise
@@ -100,11 +106,11 @@ function modal_form(sys; C1 = false)
         T = T*T2
         sysm.A .= Ab # Ab unchanged by diagonal T
     end
-    sysm, T
+    sysm, T, E
 end
 
 """
-    sysm, T = schur_form(sys)
+    sysm, T, SF = schur_form(sys)
 
 Bring `sys` to Schur form.
 
@@ -112,6 +118,7 @@ The Schur form is characterized by `A` being Schur with the real values of eigen
 ```julia
 sysm ≈ similarity_transform(sys, T)
 ```
+`SF` is the Schur-factorization of `A`.
 
 See also [`modal_form`](@ref) and [`hess_form`](@ref)
 """
@@ -124,7 +131,7 @@ function schur_form(sys)
 end
 
 """
-    sysm, T = hess_form(sys)
+    sysm, T, HF = hess_form(sys)
 
 Bring `sys` to Hessenberg form form.
 
@@ -132,6 +139,7 @@ The Hessenberg form is characterized by `A` having upper Hessenberg structure. `
 ```julia
 sysm ≈ similarity_transform(sys, T)
 ```
+`HF` is the Hessenberg-factorization of `A`.
 
 See also [`modal_form`](@ref) and [`schur_form`](@ref)
 """
diff --git a/src/descriptor.jl b/src/descriptor.jl
index 0586e2b7..8be1a062 100644
--- a/src/descriptor.jl
+++ b/src/descriptor.jl
@@ -33,6 +33,10 @@ function hinfnorm2(sys::LTISystem; kwargs...)
     DescriptorSystems.ghinfnorm(dss(ss(sys)); kwargs...)
 end
 
+function linfnorm2(sys::LTISystem; kwargs...)
+    DescriptorSystems.glinfnorm(dss(ss(sys)); kwargs...)
+end
+
 """
     n = h2norm(sys::LTISystem; kwargs...)
 
@@ -74,7 +78,7 @@ const νgap = nugap
 
 
 """
-    baltrunc2(sys::LTISystem; residual=false, n=missing, kwargs...)
+    sysr, hs = baltrunc2(sys::LTISystem; residual=false, n=missing, kwargs...)
 
 Compute the a balanced truncation of order `n` and the hankel singular values
 
@@ -87,7 +91,7 @@ function baltrunc2(sys::LTISystem; residual=false, n=missing, kwargs...)
 end
 
 """
-    baltrunc_coprime(sys; residual = false, n = missing, factorization::F = DescriptorSystems.glcf, kwargs...)
+    sysr, hs, info = baltrunc_coprime(sys; residual = false, n = missing, factorization::F = DescriptorSystems.glcf, kwargs...)
 
 Compute a balanced truncation of the left coprime factorization of `sys`.
 See [`baltrunc2`](@ref) for additional keyword-argument help.
@@ -96,15 +100,19 @@ See [`baltrunc2`](@ref) for additional keyword-argument help.
 - `factorization`: The function to perform the coprime factorization. A normalized factorization may be used by passing `RobustAndOptimalControl.DescriptorSystems.gnlcf`.
 - `kwargs`: Are passed to `DescriptorSystems.gbalmr`
 """
-function baltrunc_coprime(sys; residual=false, n=missing, factorization::F = DescriptorSystems.glcf, kwargs...) where F
-    N,M = factorization(dss(sys))
-    nu = size(N.B, 2)
-    A,E,B,C,D = DescriptorSystems.dssdata(N)
-    NM = DescriptorSystems.dss(A,E,[B M.B],C,[D M.D])
+function baltrunc_coprime(sys, info=nothing; residual=false, n=missing, factorization::F = DescriptorSystems.glcf, kwargs...) where F
+    if info !== nothing && hasproperty(info, :NM)
+        @unpack N, M, NM = info
+    else
+        N,M = factorization(dss(sys))
+        A,E,B,C,D = DescriptorSystems.dssdata(N)
+        NM = DescriptorSystems.dss(A,E,[B M.B],C,[D M.D])
+    end
     sysr, hs = DescriptorSystems.gbalmr(NM; matchdc=residual, ord=n, kwargs...)
-
+    
     A,E,B,C,D = DescriptorSystems.dssdata(DescriptorSystems.dss2ss(sysr)[1])
-
+    
+    nu = size(N.B, 2)
     BN = B[:, 1:nu]
     DN = D[:, 1:nu]
     BM = B[:, nu+1:end]
@@ -115,7 +123,7 @@ function baltrunc_coprime(sys; residual=false, n=missing, factorization::F = Des
     Br = BN  - BM * (DMi * DN)
     Dr = (DMi * DN)
 
-    ss(Ar,Br,Cr,Dr,sys.timeevol), hs
+    ss(Ar,Br,Cr,Dr,sys.timeevol), hs, (; NM, N, M)
 end
 
 
@@ -126,14 +134,18 @@ Balanced truncation for unstable models. An additive decomposition of sys into `
 
 See `baltrunc2` for other keyword arguments.
 """
-function baltrunc_unstab(sys::LTISystem; residual=false, n=missing, kwargs...)
-    stab, unstab = DescriptorSystems.gsdec(dss(sys); job="stable", kwargs...)
+function baltrunc_unstab(sys::LTISystem, info=nothing; residual=false, n=missing, kwargs...)
+    if info !== nothing && hasproperty(info, :stab)
+        @unpack stab, unstab = info
+    else
+        stab, unstab = DescriptorSystems.gsdec(dss(sys); job="stable", kwargs...)
+    end
     nx_unstab = size(unstab.A, 1)
     if n isa Integer && n < nx_unstab
         error("The model contains $(nx_unstab) poles outside the stability region, the reduced-order model must be of at least this order.")
     end
     sysr, hs = DescriptorSystems.gbalmr(stab; matchdc=residual, ord=n-nx_unstab, kwargs...)
-    ss(sysr + unstab), hs
+    ss(sysr + unstab), hs, (; stab, unstab)
 end
 
 """
diff --git a/src/glover_mcfarlane.jl b/src/glover_mcfarlane.jl
index 11023bf1..0f84421d 100644
--- a/src/glover_mcfarlane.jl
+++ b/src/glover_mcfarlane.jl
@@ -284,7 +284,7 @@ end
 
 Joint design of feedback and feedforward compensators
 ```math
-K = \\left[K_1 & K_2\\right]
+K = \\begin{bmatrix} K_1 & K_2 \\end{bmatrix}
 ```
 ```
    ┌──────┐   ┌──────┐        ┌──────┐    ┌─────┐
diff --git a/src/reduction.jl b/src/reduction.jl
index f91ee227..fc64bd5b 100644
--- a/src/reduction.jl
+++ b/src/reduction.jl
@@ -4,7 +4,7 @@ using ControlSystems: ssdata
 
 
 """
-    frequency_weighted_reduction(G, Wo, Wi; residual=true)
+    sysr, hs = frequency_weighted_reduction(G, Wo, Wi; residual=true)
 
 Find Gr such that ||Wₒ(G-Gr)Wᵢ||∞ is minimized.
 For a realtive reduction, set Wo = inv(G) and Wi = I.
@@ -113,7 +113,7 @@ function frequency_weighted_reduction(G, Wo, Wi, r=nothing; residual=true, atol=
         # determine a standard reduced system
         SV = svd!(Er)
         di2 = Diagonal(1 ./sqrt.(SV.S))
-        return ss(di2*SV.U'*Ar*SV.Vt'*di2, di2*(SV.U'*Br), (Cr*SV.Vt')*di2, Dr)
+        return ss(di2*SV.U'*Ar*SV.Vt'*di2, di2*(SV.U'*Br), (Cr*SV.Vt')*di2, Dr), Σ
     else
         L = (Y'X)\Y'
         T = X
@@ -122,7 +122,7 @@ function frequency_weighted_reduction(G, Wo, Wi, r=nothing; residual=true, atol=
         C = C*T
         D = D
     end
-    ss(A,B,C,D, G.timeevol)
+    ss(A,B,C,D, G.timeevol), Σ
 end
 
 
@@ -178,6 +178,14 @@ end
 
 @deprecate minreal2 minreal
 
+"""
+    error_bound(hs)
+
+Given a vector of Hankel singular vlues, return the theoretical error bound as a function of model order after balanced-truncation model reduction.
+(twice sum of all the removed singular values).
+"""
+error_bound(hs) = [2reverse(cumsum(reverse(hs)))[1:end-1]; 0]
+
 # slide 189 https://cscproxy.mpi-magdeburg.mpg.de/mpcsc/benner/talks/lecture-MOR.pdf
 # This implementation works, but results in a complex-valued system.
 # function model_reduction_irka(sys::AbstractStateSpace{Continuous}, r; tol = 1e-6)
diff --git a/test/test_reduction.jl b/test/test_reduction.jl
index fb4533c5..684076dc 100644
--- a/test/test_reduction.jl
+++ b/test/test_reduction.jl
@@ -9,7 +9,7 @@ sys = let
     _D = [0.0]
     ss(_A, _B, _C, _D)
 end
-sysr = frequency_weighted_reduction(sys, 1, 1, 10, residual=false)
+sysr, _ = frequency_weighted_reduction(sys, 1, 1, 10, residual=false)
 sysr2 = baltrunc(sys, n=10)[1]
 @test sysr.nx == 10
 @test hinorm(sys-sysr) < 1e-5
@@ -25,7 +25,7 @@ sys = let
     ss(_A, _B, _C, _D)
 end
 sysi = fudge_inv(sys)
-sysr = frequency_weighted_reduction(sys, sysi, 1, 5, residual=false)
+sysr, _ = frequency_weighted_reduction(sys, sysi, 1, 5, residual=false)
 sysr2 = baltrunc(sys, n=5)[1]
 @test sysr.nx == 5
 @test norm(sys-sysr) < 0.1
@@ -35,7 +35,7 @@ sysr2 = baltrunc(sys, n=5)[1]
 
 
 Wi = tf(1, [1, 1])
-sysr = frequency_weighted_reduction(sys, 1, Wi, 5, residual=false)
+sysr, _ = frequency_weighted_reduction(sys, 1, Wi, 5, residual=false)
 sysr2 = baltrunc(sys, n=5)[1]
 @test sysr.nx == 5
 @test norm(sys-sysr) < 0.1
@@ -47,7 +47,7 @@ sysr2 = baltrunc(sys, n=5)[1]
 
 Wo = sysi
 Wi = tf(1, [1, 1])
-sysr = frequency_weighted_reduction(sys, Wo, Wi, 5, residual=false)
+sysr, _ = frequency_weighted_reduction(sys, Wo, Wi, 5, residual=false)
 sysr2 = baltrunc(sys, n=5)[1]
 @test sysr.nx == 5
 @test_broken norm(sys-sysr) < 0.1
@@ -58,7 +58,7 @@ sysr2 = baltrunc(sys, n=5)[1]
 
 # residual
 sysi = fudge_inv(sys)
-sysr = frequency_weighted_reduction(sys, sysi, 1, 3, residual=true)
+sysr, _ = frequency_weighted_reduction(sys, sysi, 1, 3, residual=true)
 sysr2 = baltrunc(sys, n=3, residual=true)[1]
 @test sysr.nx == 3
 @test norm(sys-sysr, Inf) < 3
@@ -72,7 +72,7 @@ sysr2 = baltrunc(sys, n=3, residual=true)[1]
 G = tf([8, 6, 2], [1, 4, 5, 2])
 Wi = tf(1, [1, 3])
 Wo = tf(1, [1, 4])
-sysr = frequency_weighted_reduction(ss(G), Wo, Wi, 1, residual=true)
+sysr, _ = frequency_weighted_reduction(ss(G), Wo, Wi, 1, residual=true)
 @test isstable(sysr)
 
 @test tf(sysr) ≈ tf([2.398, 1.739], [1, 1.739]) rtol=1e-1
@@ -100,7 +100,7 @@ Wo = Wi = ss(Aw,Dw,Cw,Dw)
 
 
 errors = map(1:3) do r
-    Gr = frequency_weighted_reduction(G, Wo, Wi, r, residual=false)
+    Gr, _ = frequency_weighted_reduction(G, Wo, Wi, r, residual=false)
     hinfnorm2(Wo*(G-Gr)*Wi)[1]
 end
 
@@ -110,7 +110,7 @@ end
 
 
 errors = map(1:3) do r
-    Gr = frequency_weighted_reduction(G, Wo, Wi, r, residual=true)
+    Gr, _ = frequency_weighted_reduction(G, Wo, Wi, r, residual=true)
     hinfnorm2(Wo*(G-Gr)*Wi)[1]
 end