diff --git a/libsolutil/LP.cpp b/libsolutil/LP.cpp
index d0cf68b6f..292912726 100644
--- a/libsolutil/LP.cpp
+++ b/libsolutil/LP.cpp
@@ -273,6 +273,13 @@ void LPSolver::combineSubProblems(size_t _combineInto, size_t _combineFrom)
 	SubProblem const& combineFrom = *m_subProblems[_combineFrom];
 
 	size_t varShift = combineInto.variables.size();
+#ifdef SPARSE
+	size_t rowShift = combineInto.factors.rows();
+
+	for (size_t row = 0; row < combineFrom.factors.rows(); row++)
+		for (auto&& [col, v]: combineFrom.factors.iterateRow(row))
+			combineinto.factors.insert(row + rowShift, col + colShift, move(v));
+#else
 	size_t rowShift = combineInto.factors.size();
 	size_t newRowLength = combineInto.variables.size() + combineFrom.variables.size();
 	for (LinearExpression& row: combineInto.factors)
@@ -285,6 +292,7 @@ void LPSolver::combineSubProblems(size_t _combineInto, size_t _combineFrom)
 			shiftedRow[varShift + index] = f;
 		combineInto.factors.emplace_back(move(shiftedRow));
 	}
+#endif
 	combineInto.variables += combineFrom.variables;
 	for (auto const& index: combineFrom.variablesPotentiallyOutOfBounds)
 		combineInto.variablesPotentiallyOutOfBounds.insert(index + varShift);
diff --git a/libsolutil/LP.h b/libsolutil/LP.h
index c06c6b1d4..df5ba7da2 100644
--- a/libsolutil/LP.h
+++ b/libsolutil/LP.h
@@ -17,6 +17,9 @@
 // SPDX-License-Identifier: GPL-3.0
 #pragma once
 
+// use sparse matrices
+#define SPARSE 1
+
 #include <libsolutil/Numeric.h>
 #include <libsolutil/LinearExpression.h>
 
@@ -205,7 +208,11 @@ private:
 		/// Set to true on "check". Needs a copy for adding a constraint or bound if set to true.
 		bool sealed = false;
 		std::optional<LPResult> result = std::nullopt;
+#ifdef SPARSE
+		SparseMatrix factors;
+#else
 		std::vector<LinearExpression> factors;
+#endif
 		std::vector<Variable> variables;
 		std::set<size_t> variablesPotentiallyOutOfBounds;
 		/// Variable index to constraint it controls.
diff --git a/libsolutil/LinearExpression.cpp b/libsolutil/LinearExpression.cpp
index b89e8038d..b8c61fcc5 100644
--- a/libsolutil/LinearExpression.cpp
+++ b/libsolutil/LinearExpression.cpp
@@ -21,6 +21,37 @@
 using namespace solidity::util;
 using namespace std;
 
+SparseMatrix::SparseMatrixIterator SparseMatrix::IteratorCombiner::begin()
+{
+	return SparseMatrixIterator(
+		m_isRow ? m_matrix.m_row_start[m_RowOrColumn] : m_matrix.m_col_start[m_RowOrColumn],
+		m_isRow
+	);
+}
+
+SparseMatrix::SparseMatrixIterator SparseMatrix::IteratorCombiner::end()
+{
+	return SparseMatrixIterator(nullptr, m_isRow);
+}
+
+SparseMatrix::IteratorCombiner SparseMatrix::enumerateColumn(size_t _column)
+{
+	return IteratorCombiner{
+		_column,
+		false,
+		*this
+	};
+}
+
+SparseMatrix::IteratorCombiner SparseMatrix::enumerateRow(size_t _row)
+{
+	return IteratorCombiner{
+		_row,
+		true,
+		*this
+	};
+}
+
 void SparseMatrix::multiplyRowByFactor(size_t _row, rational const& _factor)
 {
 	Entry* e = m_row_start[_row];
@@ -68,7 +99,7 @@ void SparseMatrix::appendRow(LinearExpression const& _entries)
 	for (auto&& [i, v]: _entries.enumerate()) {
 		if (!v)
 			continue;
-		appendToRow(row_nr, i, move(v));
+		prependInRow(nullptr, row_nr, i, move(v));
 	}
 }
 
diff --git a/libsolutil/LinearExpression.h b/libsolutil/LinearExpression.h
index 65586555d..9c57f0c25 100644
--- a/libsolutil/LinearExpression.h
+++ b/libsolutil/LinearExpression.h
@@ -1,4 +1,4 @@
-/*
+/*
 	This file is part of solidity.
 
 	solidity is free software: you can redistribute it and/or modify
@@ -183,6 +183,7 @@ private:
 	std::vector<rational> factors;
 };
 
+
 class SparseMatrix
 {
 public:
@@ -198,19 +199,64 @@ public:
 		Entry* prev_in_col;
 		Entry* next_in_col;
 	};
+	struct SparseMatrixIterator
+	{
+		using iterator_category = std::forward_iterator_tag;
+		using difference_type   = std::ptrdiff_t;
+		using value_type        = Entry;
+		using pointer           = Entry*;
+		using reference         = Entry&;
+
+		SparseMatrixIterator(pointer _ptr, bool _isRow): m_ptr(_ptr), m_isRow(_isRow) {}
+
+		reference operator*() const { return *m_ptr; }
+		pointer operator->() { return m_ptr; }
+		SparseMatrixIterator& operator++()
+		{
+			m_ptr = m_isRow ? m_ptr->next_in_row : m_pt->next_in_col;
+			return *this;
+		}
+		SparseMatrixIterator operator++(int) { SparseMatrixIterator tmp = *this; ++(*this); return tmp; }
+		friend bool operator==(SparseMatrixIterator const& _a, SparseMatrixIterator const& _b)
+		{
+			return _a.m_ptr == _b.m_ptr && _a.m_isRow == _b.m_isRow;
+		}
+		friend bool operator!=(SparseMatrixIterator const& _a, SparseMatrixIterator const& _b)
+		{
+			return _a.m_ptr != _b.m_ptr || _a.m_isRow != _b.m_isRow;
+		}
+
+	private:
+		Entry* m_ptr;
+		bool m_isRow;
+	};
+	struct IteratorCombiner
+	{
+		size_t m_RowOrColumn;
+		bool m_isRow;
+		SparseMatrix& m_matrix;
+		SparseMatrixIterator begin();
+		SparseMatrixIterator end();
+	};
+
+	size_t rows() const { return m_row_start.size(); }
+	size_t columns() const { return m_col_start.size(); }
+
 	/// @returns (i, v) for all non-zero v in the column _column
-	void enumerateColumn(size_t _column) const;
+	void enumerateColumn(size_t _column);
 	/// @returns (i, v) for all non-zero v in the row _row
-	void enumerateRow(size_t _row) const;
+	void enumerateRow(size_t _row);
 	void multiplyRowByFactor(size_t _row, rational const& _factor);
 	void addMultipleOfRow(size_t _sourceRow, size_t _targetRow, rational const& _factor);
 	rational entry(size_t _row, size_t _column) const;
+	void insert(size_t _row, size_t _column, rational _value);
 
 	void appendRow(LinearExpression const& _entries);
 
 private:
 
 	void remove(Entry& _entry);
+	/// Prepends a new entry before the given element or at end of row if nullptr.
 	Entry* prependInRow(Entry* _successor, size_t _row, size_t _column, rational _value);
 	void adjustColumnProperties(Entry& _entry);