diff --git a/cpp/Matrix.cpp b/cpp/Matrix.cpp
index e7f87098d..04ab4b2ff 100644
--- a/cpp/Matrix.cpp
+++ b/cpp/Matrix.cpp
@@ -150,16 +150,18 @@ Vector Vector_(const Matrix& A)
 }
 
 /* ************************************************************************* */
-Vector column(const Matrix& A, size_t j) {
-	if (j>=A.size2())
-		throw invalid_argument("Column index out of bounds!");
+Vector column_(const Matrix& A, size_t j) {
+//	if (j>=A.size2())
+//		throw invalid_argument("Column index out of bounds!");
+
+	return column(A,j); // real boost version
 
 	// TODO: improve this
-	size_t m = A.size1();
-	Vector a(m);
-	for (size_t i=0; i<m; ++i)
-		a(i) = A(i,j);
-	return a;
+//	size_t m = A.size1();
+//	Vector a(m);
+//	for (size_t i=0; i<m; ++i)
+//		a(i) = A(i,j);
+//	return a;
 }
 
 /* ************************************************************************* */
@@ -336,7 +338,7 @@ weighted_eliminate(Matrix& A, Vector& b, const Vector& sigmas) {
 	// Then update A and b by substituting x with d-rS, zero-ing out x's column.
 	for (size_t j=0; j<n; ++j) {
 		// extract the first column of A
-		Vector a(column(A, j)); // ublas::matrix_column is slower !
+		Vector a(column_(A, j)); // ublas::matrix_column is slower !
 
 		// Calculate weighted pseudo-inverse and corresponding precision
 		double precision = weightedPseudoinverse(a, weights, pseudo);
@@ -347,7 +349,7 @@ weighted_eliminate(Matrix& A, Vector& b, const Vector& sigmas) {
 		// create solution and copy into r
 		Vector r(basis(n, j));
 		for (size_t j2=j+1; j2<n; ++j2)
-			r(j2) = inner_prod(pseudo, ublas::matrix_column<Matrix>(A, j2));
+			r(j2) = inner_prod(pseudo, ublas::matrix_column<Matrix>(A, j2)); // TODO: don't use ublas
 
 		// create the rhs
 		double d = inner_prod(pseudo, b);
@@ -359,7 +361,7 @@ weighted_eliminate(Matrix& A, Vector& b, const Vector& sigmas) {
 		// exit after rank exhausted
 		if (results.size()>=maxRank) break;
 
-		// update A, b, expensive, suing outer product
+		// update A, b, expensive, using outer product
 		// A' \define A_{S}-a*r and b'\define b-d*a
 		updateAb(A, b, j, a, r, d);
 	}
diff --git a/cpp/Matrix.h b/cpp/Matrix.h
index 723523013..813e15d8e 100644
--- a/cpp/Matrix.h
+++ b/cpp/Matrix.h
@@ -144,11 +144,12 @@ Matrix sub(const Matrix& A, size_t i1, size_t i2, size_t j1, size_t j2);
 
 /**
  * extracts a column from a matrix
+ * NOTE: using this without the underscore is the ublas version!
  * @param matrix to extract column from
  * @param index of the column
  * @return the column in vector form
  */
-Vector column(const Matrix& A, size_t j);
+Vector column_(const Matrix& A, size_t j);
 
 /**
  * extracts a row from a matrix
diff --git a/cpp/NoiseModel.cpp b/cpp/NoiseModel.cpp
index bd50756c9..f9c985a99 100644
--- a/cpp/NoiseModel.cpp
+++ b/cpp/NoiseModel.cpp
@@ -188,7 +188,10 @@ namespace gtsam {
 			// Then update A and b by substituting x with d-rS, zero-ing out x's column.
 			for (size_t j=0; j<n; ++j) {
 				// extract the first column of A
-				Vector a(column(Ab, j)); // ublas::matrix_column is slower ! TODO Really, why ????
+				// ublas::matrix_column is slower ! TODO Really, why ????
+				//  AGC: if you use column() you will automatically call ublas, use
+				//      column_() to actually use the one in our library
+				Vector a(column(Ab, j));
 
 				// Calculate weighted pseudo-inverse and corresponding precision
 				double precision = weightedPseudoinverse(a, weights, pseudo);
diff --git a/cpp/testMatrix.cpp b/cpp/testMatrix.cpp
index 0b26f75d4..894d974dc 100644
--- a/cpp/testMatrix.cpp
+++ b/cpp/testMatrix.cpp
@@ -153,15 +153,15 @@ TEST( matrix, column )
 		0., -1.,  0.,  1.,  0.,  0.,  0.3,
 		1.,  0.,  0.,  0., -1.,  0.,  0.2,
 		0.,  1.,  0.,  0.,  0., -1., -0.1);
-	Vector a1 = column(A, 0);
+	Vector a1 = column_(A, 0);
 	Vector exp1 = Vector_(4, -1., 0., 1., 0.);
 	CHECK(assert_equal(a1, exp1));
 
-	Vector a2 = column(A, 3);
+	Vector a2 = column_(A, 3);
 	Vector exp2 = Vector_(4,  0., 1., 0., 0.);
 	CHECK(assert_equal(a2, exp2));
 
-	Vector a3 = column(A, 6);
+	Vector a3 = column_(A, 6);
 	Vector exp3 = Vector_(4, -0.2, 0.3, 0.2, -0.1);
 	CHECK(assert_equal(a3, exp3));
 }
diff --git a/cpp/timeMatrix.cpp b/cpp/timeMatrix.cpp
index 38ebc9579..afcd64996 100644
--- a/cpp/timeMatrix.cpp
+++ b/cpp/timeMatrix.cpp
@@ -5,11 +5,13 @@
  */
 
 #include <iostream>
+#include <boost/numeric/ublas/matrix_proxy.hpp>
 #include <boost/timer.hpp>
 #include "Matrix.h"
 
 using namespace std;
 using namespace gtsam;
+namespace ublas = boost::numeric::ublas;
 
 /*
  * Results:
@@ -114,13 +116,41 @@ double timeVScaleRow(size_t m, size_t n, size_t reps) {
 	{
 		boost::timer t;
 		for (int i=0; i<reps; ++i)
-			Matrix result = vector_scale(V,M);
+			result = vector_scale(V,M);
 		elapsed = t.elapsed();
 	}
 
 	return elapsed;
 }
 
+/**
+ * Results:
+ * Alex's Machine
+ *  - ublas matrix_column  : 4.63 sec
+ *  - naive implementation : 4.70 sec
+ */
+double timeColumn(size_t reps) {
+	// create a matrix
+	size_t m = 100; size_t n = 100;
+	Matrix M(m, n);
+	for (int i=0; i<m; ++i)
+			for (int j=0; j<n; ++j)
+				M(i,j) = 2*i+j;
+
+	// extract a column
+	double elapsed;
+	Vector result;
+	{
+		boost::timer t;
+		for (size_t i=0; i<reps; ++i)
+			for (size_t j = 0; j<n; ++j)
+				result = ublas::matrix_column<Matrix>(M, j);
+				//result = column(M, j);
+		elapsed = t.elapsed();
+	}
+	return elapsed;
+}
+
 int main(int argc, char ** argv) {
 
 	// Time collect()
@@ -143,5 +173,11 @@ int main(int argc, char ** argv) {
 	double vsRow_time = timeVScaleRow(m1, n1, reps1);
 	cout << "Elapsed time for vector_scale(row)    [(" << m1 << ", " << n1 << ") matrix] : " << vsRow_time << endl;
 
+	// Time column_() NOTE: using the gtsam version
+	cout << "Starting column_() Timing" << endl;
+	size_t reps2 = 200000;
+	double column_time = timeColumn(reps2);
+	cout << "Time: " << column_time << " sec" << endl;
+
 	return 0;
 }