diff --git a/src/ReservoirComputing.jl b/src/ReservoirComputing.jl
index 153ddd72..a85e2015 100644
--- a/src/ReservoirComputing.jl
+++ b/src/ReservoirComputing.jl
@@ -3,7 +3,7 @@ module ReservoirComputing
 using Adapt: adapt
 using CellularAutomata: CellularAutomaton
 using Compat: @compat
-using LinearAlgebra: eigvals, mul!, I
+using LinearAlgebra: eigvals, mul!, I, qr
 using NNlib: fast_act, sigmoid
 using Random: Random, AbstractRNG
 using Reexport: Reexport, @reexport
diff --git a/src/train/linear_regression.jl b/src/train/linear_regression.jl
index c66ff104..429b816f 100644
--- a/src/train/linear_regression.jl
+++ b/src/train/linear_regression.jl
@@ -30,21 +30,11 @@ function StandardRidge()
 end
 
 function train(sr::StandardRidge, states::AbstractArray, target_data::AbstractArray)
-    #A = states * states' + sr.reg * I
-    #b = states * target_data
-    #output_layer = (A \ b)'
-
-    if size(states, 2) != size(target_data, 2)
-        throw(DimensionMismatch("\n" *
-                                "\n" *
-                                "  - Number of columns in `states`: $(size(states, 2))\n" *
-                                "  - Number of columns in `target_data`: $(size(target_data, 2))\n" *
-                                "The dimensions of `states` and `target_data` must align for training." *
-                                "\n"
-        ))
-    end
-
-    output_layer = Matrix(((states * states' + sr.reg * I) \
-                           (states * target_data'))')
+    n_states = size(states, 1)
+    A = [states'; sqrt(sr.reg) * I(n_states)]
+    b = [target_data'; zeros(n_states, size(target_data, 1))]
+    F = qr(A)
+    Wt = F \ b
+    output_layer = Matrix(Wt')
     return OutputLayer(sr, output_layer, size(target_data, 1), target_data[:, end])
 end