diff --git a/src/Orbit-wise.jl b/src/Orbit-wise.jl
index a305db4..bf85296 100644
--- a/src/Orbit-wise.jl
+++ b/src/Orbit-wise.jl
@@ -53,13 +53,18 @@ end
 
 include("OrbitDecomposition.jl")
 
-function sparsify{T}(U::AbstractArray{T}, eps=eps(T), check=true)
+dens(M::SparseMatrixCSC) = length(M.nzval)/length(M)
+dens(M::AbstractArray) = sum(abs.(M) .!= 0)/length(M)
+
+function sparsify{T}(U::AbstractArray{T}, check=true)
     W = deepcopy(U)
-    W[abs.(W) .< eps] = zero(T)
+    W[abs.(W) .< eps(T)] .= zero(T)
     if check && rank(W) != rank(U)
-        warn("Sparsification would decrease the rank!")
+        info("Sparsification would decrease the rank!")
         W = U
-    end
+    else
+      info("Sparsified:", rpad(dens(U), 10), "\t", rpad(dens(W),10))
+   end
     W = sparse(W)
     dropzeros!(W)
     return W
@@ -78,15 +83,15 @@ function init_orbit_data(logger, sett::Settings; radius=2)
    return 0
 end
 
-function transform(U::AbstractArray, V::AbstractArray; sparse=false)
-   w = U'*V*U
+function transform(U::AbstractArray, V::AbstractArray; sparse=true)
    if sparse
-      w = sparsify(w)
+      return sparsify(U'*V*U)
+   else
+      return U'*V*U
    end
-   return w
 end
 
-A(data::OrbitData, π, t) = data.dims[π]*transform(data.Us[π], data.cnstr[t], sparse=true)
+A(data::OrbitData, π, t) = data.dims[π]*transform(data.Us[π], data.cnstr[t])
 
 function constrLHS(m::JuMP.Model, data::OrbitData, t)
     l = endof(data.Us)
diff --git a/src/OrbitDecomposition.jl b/src/OrbitDecomposition.jl
index b8e7eb7..30b7ec3 100644
--- a/src/OrbitDecomposition.jl
+++ b/src/OrbitDecomposition.jl
@@ -123,7 +123,7 @@ function reconstruct_sol{T<:GroupElem, S<:AbstractArray}(mreps::Dict{T, S},
     for g in keys(mreps)
         A, B = mreps[g], mreps[inv(g)]
         for π in 1:length(Us)
-            recP .+= dims[π].* (A * Us[π]*Ps[π]*Ust[π] * B)
+            recP .+= sparsify(dims[π].* (A * Us[π]*Ps[π]*Ust[π] * B))
         end
     end
     recP .*= 1/length(keys(mreps))
@@ -140,19 +140,6 @@ function Cstar_repr(x::GroupRingElem, mreps::Dict)
     return res
 end
 
-dens(M::SparseMatrixCSC) = length(M.nzval)/length(M)
-dens(M::AbstractArray) = sum(abs.(M) .!= 0)/length(M)
-
-function sparsify2(M::AbstractArray)
-   println("Density before sparsification: \t$(dens(M))")
-   M = deepcopy(M)
-   M[M .< eps(eltype(M))] .= zero(eltype(M))
-   M = sparse(M)
-   dropzeros!(M)
-   println("Density after sparsification: \t$(dens(M))")
-   return M
-end
-
 function orthSVD(M::AbstractMatrix)
     M = full(M)
     fact = svdfact(M)
@@ -205,6 +192,6 @@ function compute_orbit_data{T<:GroupElem}(logger, name::String, G::Nemo.Group, S
    info(logger, "dimensions = $dimensions")
    @assert dot(multiplicities, dimensions) == sizes[radius]
 
-   save(joinpath(name, "U_pis.jld"), "Uπs", Uπs, "spUπs", sparsify2.(Uπs), "dims", dimensions)
+   save(joinpath(name, "U_pis.jld"), "Uπs", Uπs, "spUπs", sparsify.(Uπs), "dims", dimensions)
    return 0
 end