diff --git a/lib/NumbersPHP/Matrix.php b/lib/NumbersPHP/Matrix.php
index 29875b8..ddd6588 100644
--- a/lib/NumbersPHP/Matrix.php
+++ b/lib/NumbersPHP/Matrix.php
@@ -192,6 +192,74 @@ public static function determinant($matrix)
         }
         return $determinant;
     }
+    
+    /**
+    * Inverts a given matrix.
+    * 
+    * For a given matrix A, inverse B is given as AB = BA = I.
+    *
+    * @param array $A matrix to invert
+    * @return $Inv inverted matrix
+    * @throws \Exception 
+    */
+    function invert($A)
+    {
+	    if (!self::isSquare($matrix)) {
+            throw new \Exception('Matrix must be square');
+        }
+        
+	    $n = count($A);
+	    // get and append identity matrix
+	    $I = identity_matrix($n);
+	    for ($i = 0; $i < $n; ++ $i) {
+    		$A[$i] = array_merge($A[$i], $I[$i]);
+    	}
+    	
+	    // forward run
+	    for ($j = 0; $j < $n-1; ++ $j) {
+    		// for all remaining rows (diagonally)
+		    for ($i = $j+1; $i < $n; ++ $i) {
+    			// if the value is not already 0
+			    if ($A[$i][$j] !== 0) {
+    				// adjust scale to pivot row
+				    // subtract pivot row from current
+				    $scalar = $A[$j][$j] / $A[$i][$j];
+				    for ($jj = $j; $jj < $n*2; ++ $jj) {
+    					$A[$i][$jj] *= $scalar;
+					    $A[$i][$jj] -= $A[$j][$jj];
+				    }
+			    }
+		    }
+	    }
+	    // reverse run
+	    for ($j = $n-1; $j > 0; -- $j) {
+    		for ($i = $j-1; $i >= 0; -- $i) {
+			    if ($A[$i][$j] !== 0) {
+    				$scalar = $A[$j][$j] / $A[$i][$j];
+				    for ($jj = $i; $jj < $n*2; ++ $jj) {
+    					$A[$i][$jj] *= $scalar;
+					    $A[$i][$jj] -= $A[$j][$jj];
+				    }
+			    }
+		    }
+	    }
+	    // last run to make all diagonal 1s
+	    /// @note this can be done in last iteration (i.e. reverse run) too!
+	    for ($j = 0; $j < $n; ++ $j) {
+    		if ($A[$j][$j] !== 1) {
+			    $scalar = 1 / $A[$j][$j];
+			    for ($jj = $j; $jj < $n*2; ++ $jj) {
+    				$A[$j][$jj] *= $scalar;
+			    }
+		    }
+	    }
+	    // take out the matrix inverse to return
+	    $Inv = array();
+	    for ($i = 0; $i < $n; ++ $i) {
+    		$Inv[$i] = array_slice($A[$i], $n);
+    	}
+    	return $Inv;
+    }
 
     /**
      * Returns a LUP decomposition of the given matrix such that: