Fix numerical stability issues in SVD implementation (#73)

* Fix numerical stability issues in SVD implementation

- Add EPSILON constant for safe numerical comparisons
- Guard Vector#normalize against zero/near-zero magnitude
- Guard SVD angle calculation against division by zero
- Guard SVD singular value sqrt against negative values
- Guard content_node weighted_total against division by zero
- Add nil checks in LSI build_index for edge cases
- Add tests for zero vectors, near-singular matrices, edge cases

Fixes #63

* Fix RuboCop offenses in SVD implementation

* Address code review feedback

- Move EPSILON constant to Vector::EPSILON to avoid polluting global namespace
- Use Vector::EPSILON in content_node.rb for consistency
- Remove unrelated .claude/settings.local.json from PR
- Add .claude/ to .gitignore

* Address additional code review feedback

- Change normalize guard from `< EPSILON` to `<= 0.0` - only return zero
  vector for actual zero magnitude, not near-zero values
- Use `val.abs` with EPSILON threshold for singular values to handle
  small negative values from floating point errors more robustly
- Update test to reflect correct behavior: near-zero vectors normalize
  to unit vectors, only true zero vectors return zero

* Fix RuboCop: add empty line before assertions

* Address third round of code review feedback

- Rename test to test_normalize_near_zero_vector_normalizes_correctly
- Simplify singular value sqrt to [val, 0.0].max for uniform handling
- Preserve sign relationship in content_node divisor fallback

* Remove redundant comments

* Use Vector::EPSILON in test for maintainability

* Flatten nested conditional in content_node divisor logic

Author: cardmagic

.gitignore

@@ -13,3 +13,6 @@ pkg/
 
 # OS files
 .DS_Store
+
+# Claude Code local settings
+.claude/

lib/classifier/extensions/vector.rb

@@ -14,7 +14,9 @@ def sum_with_identity(identity = 0.0, &)
   end
 end
 
-module VectorExtensions
+class Vector
+  EPSILON = 1e-10
+
   def magnitude
     sum_of_squares = 0.to_r
     size.times do |i|
@@ -24,19 +26,17 @@ def magnitude
   end
 
   def normalize
+    magnitude_value = magnitude
+    return Vector[*Array.new(size, 0.0)] if magnitude_value <= 0.0
+
     normalized_values = []
-    magnitude_value = magnitude.to_r
     size.times do |i|
       normalized_values << (self[i] / magnitude_value)
     end
     Vector[*normalized_values]
   end
 end
 
-class Vector
-  include VectorExtensions
-end
-
 class Matrix
   def self.diag(diagonal_elements)
     Matrix.diagonal(*diagonal_elements)
@@ -62,10 +62,15 @@ def SV_decomp(max_sweeps = 20)
         (1..(q_rotation_matrix.row_size - 1)).each do |col|
           next if row == col
 
-          angle = Math.atan((2.to_r * q_rotation_matrix[row,
-                                                        col]) / (q_rotation_matrix[row,
-                                                                                   row] - q_rotation_matrix[col,
-                                                                                                            col])) / 2.0
+          numerator = 2.0 * q_rotation_matrix[row, col]
+          denominator = q_rotation_matrix[row, row] - q_rotation_matrix[col, col]
+
+          angle = if denominator.abs < Vector::EPSILON
+                    numerator >= 0 ? Math::PI / 4.0 : -Math::PI / 4.0
+                  else
+                    Math.atan(numerator / denominator) / 2.0
+                  end
+
           cosine = Math.cos(angle)
           sine = Math.sin(angle)
           rotation_matrix = Matrix.identity(q_rotation_matrix.row_size)
@@ -89,11 +94,12 @@ def SV_decomp(max_sweeps = 20)
       break if (sum_of_differences <= 0.001 && iteration_count > 1) || iteration_count >= max_sweeps
     end
 
-    singular_values = []
-    q_rotation_matrix.row_size.times do |r|
-      singular_values << Math.sqrt(q_rotation_matrix[r, r].to_f)
+    singular_values = q_rotation_matrix.row_size.times.map do |r|
+      Math.sqrt([q_rotation_matrix[r, r].to_f, 0.0].max)
     end
-    u_matrix = (row_size >= column_size ? self : trans) * v_matrix * Matrix.diagonal(*singular_values).inverse
+
+    safe_singular_values = singular_values.map { |v| [v, Vector::EPSILON].max }
+    u_matrix = (row_size >= column_size ? self : trans) * v_matrix * Matrix.diagonal(*safe_singular_values).inverse
     [u_matrix, v_matrix, singular_values]
   end
 

lib/classifier/lsi.rb

@@ -143,8 +143,13 @@ def build_index(cutoff = 0.75)
         ntdm = build_reduced_matrix(tdm, cutoff)
 
         ntdm.row_size.times do |col|
-          doc_list[col].lsi_vector = ntdm.column(col) if doc_list[col]
-          doc_list[col].lsi_norm = ntdm.column(col).normalize if doc_list[col]
+          next unless doc_list[col]
+
+          column = ntdm.column(col)
+          next unless column
+
+          doc_list[col].lsi_vector = column
+          doc_list[col].lsi_norm = column.normalize
         end
       end
 

lib/classifier/lsi/content_node.rb

@@ -60,7 +60,10 @@ def raw_vector_with(word_list)
           val = term_over_total * Math.log(term_over_total)
           weighted_total += val unless val.nan?
         end
-        vec = vec.collect { |val| Math.log(val + 1) / -weighted_total }
+
+        sign = weighted_total.negative? ? 1.0 : -1.0
+        divisor = sign * [weighted_total.abs, Vector::EPSILON].max
+        vec = vec.collect { |val| Math.log(val + 1) / divisor }
       end
 
       if Classifier::LSI.gsl_available

test/extensions/word_hash_test.rb

@@ -115,6 +115,30 @@ def test_normalize_preserves_direction
     assert_in_delta 1.0, normalized[0], 0.001
     assert_in_delta 0.0, normalized[1], 0.001
   end
+
+  def test_normalize_zero_vector_returns_zero_vector
+    vec = Vector[0, 0, 0]
+    normalized = vec.normalize
+
+    assert_in_delta 0.0, normalized[0], 0.001
+    assert_in_delta 0.0, normalized[1], 0.001
+    assert_in_delta 0.0, normalized[2], 0.001
+  end
+
+  def test_normalize_near_zero_vector_normalizes_correctly
+    # Near-zero vectors should still normalize to unit vectors
+    # Only actual zero vectors return zero
+    small = Vector::EPSILON * 10
+    vec = Vector[small, small, small]
+    normalized = vec.normalize
+
+    # Should normalize to [1/sqrt(3), 1/sqrt(3), 1/sqrt(3)]
+    expected = 1.0 / Math.sqrt(3)
+
+    assert_in_delta expected, normalized[0], 0.001
+    assert_in_delta expected, normalized[1], 0.001
+    assert_in_delta expected, normalized[2], 0.001
+  end
 end
 
 class MatrixExtensionsTest < Minitest::Test
@@ -137,4 +161,31 @@ def test_matrix_element_assignment
 
     assert_equal 99, matrix[0, 1]
   end
+
+  def test_svd_basic
+    matrix = Matrix[[1, 0], [0, 1], [0, 0]]
+    _u, _v, s = matrix.SV_decomp
+
+    assert_equal 2, s.size
+    assert(s.all? { |val| val >= 0 })
+  end
+
+  def test_svd_with_zero_rows
+    # Matrix with linearly dependent rows that could cause zero singular values
+    matrix = Matrix[[1, 1], [1, 1], [0, 0]]
+    _u, _v, s = matrix.SV_decomp
+
+    # Should not raise an error
+    assert_equal 2, s.size
+  end
+
+  def test_svd_handles_near_singular_matrix
+    # Near-singular matrix that previously caused Math::DomainError
+    matrix = Matrix[[1e-10, 0], [0, 1e-10]]
+
+    # Should not raise an error
+    _u, _v, s = matrix.SV_decomp
+
+    assert_equal 2, s.size
+  end
 end

test/lsi/lsi_test.rb

@@ -263,4 +263,48 @@ def test_categories_for_nonexistent_item
 
     assert_empty result, 'Should return empty array for nonexistent item'
   end
+
+  # Numerical stability tests
+
+  def test_identical_documents
+    lsi = Classifier::LSI.new auto_rebuild: false
+
+    # Identical documents could cause zero singular values
+    lsi.add_item 'Exactly the same text', 'A'
+    lsi.add_item 'Exactly the same text', 'A'
+    lsi.add_item 'Different content here', 'B'
+
+    # Should handle gracefully without crashing
+    lsi.build_index
+
+    refute_predicate lsi, :needs_rebuild?
+  end
+
+  def test_single_word_documents
+    lsi = Classifier::LSI.new auto_rebuild: false
+
+    # Very short documents could cause edge cases
+    lsi.add_item 'dog', 'Animal'
+    lsi.add_item 'cat', 'Animal'
+    lsi.add_item 'car', 'Vehicle'
+
+    # Should handle gracefully
+    lsi.build_index
+
+    refute_predicate lsi, :needs_rebuild?
+  end
+
+  def test_empty_word_hash_handling
+    lsi = Classifier::LSI.new auto_rebuild: false
+
+    # Documents with only stop words result in empty word hashes
+    lsi.add_item 'the a an', 'StopWords'
+    lsi.add_item 'Dogs are great', 'Animal'
+    lsi.add_item 'Cats are nice', 'Animal'
+
+    # Should handle gracefully
+    lsi.build_index
+
+    refute_predicate lsi, :needs_rebuild?
+  end
 end