Implement left and right null

Author: mtfishman

src/KroneckerArrays.jl

@@ -381,13 +381,16 @@ using MatrixAlgebraKit:
   TruncationStrategy,
   default_eig_algorithm,
   default_eigh_algorithm,
+  default_qr_algorithm,
   eig_full!,
   eig_trunc!,
   eig_vals!,
   eigh_full!,
   eigh_trunc!,
   eigh_vals!,
   initialize_output,
+  left_null!,
+  right_null!,
   truncate!
 
 struct KroneckerAlgorithm{A,B} <: AbstractAlgorithm
@@ -461,4 +464,33 @@ function MatrixAlgebraKit.eigh_vals!(a::KroneckerMatrix, F, alg::KroneckerAlgori
   return F
 end
 
+function MatrixAlgebraKit.default_qr_algorithm(a::KroneckerMatrix; kwargs...)
+  return KroneckerAlgorithm(
+    default_qr_algorithm(a.a; kwargs...), default_qr_algorithm(a.b; kwargs...)
+  )
+end
+function MatrixAlgebraKit.default_lq_algorithm(a::KroneckerMatrix; kwargs...)
+  return KroneckerAlgorithm(
+    default_lq_algorithm(a.a; kwargs...), default_lq_algorithm(a.b; kwargs...)
+  )
+end
+
+function MatrixAlgebraKit.initialize_output(f::typeof(left_null!), a::KroneckerMatrix)
+  return initialize_output(f, a.a) ⊗ initialize_output(f, a.b)
+end
+function MatrixAlgebraKit.left_null!(a::KroneckerMatrix, F; kwargs...)
+  left_null!(a.a, F.a; kwargs...)
+  left_null!(a.b, F.b; kwargs...)
+  return F
+end
+
+function MatrixAlgebraKit.initialize_output(f::typeof(right_null!), a::KroneckerMatrix)
+  return initialize_output(f, a.a) ⊗ initialize_output(f, a.b)
+end
+function MatrixAlgebraKit.right_null!(a::KroneckerMatrix, F; kwargs...)
+  right_null!(a.a, F.a; kwargs...)
+  right_null!(a.b, F.b; kwargs...)
+  return F
+end
+
 end

test/test_matrixalgebrakit.jl

@@ -1,5 +1,5 @@
 using KroneckerArrays: ⊗
-using LinearAlgebra: Hermitian, diag
+using LinearAlgebra: Hermitian, diag, norm
 using MatrixAlgebraKit:
   eig_full,
   eig_trunc,
@@ -24,24 +24,35 @@ using MatrixAlgebraKit:
 using Test: @test, @test_throws, @testset
 
 @testset "MatrixAlgebraKit" begin
-  x = randn(2, 2)
-  y = randn(3, 3)
-  a = x ⊗ y
-  ah = Hermitian(x) ⊗ Hermitian(y)
+  elt = Float32
 
+  a = randn(elt, 2, 2) ⊗ randn(elt, 3, 3)
   d, v = eig_full(a)
   @test a * v ≈ v * d
 
+  a = randn(elt, 2, 2) ⊗ randn(elt, 3, 3)
   @test_throws MethodError eig_trunc(a)
 
+  a = randn(elt, 2, 2) ⊗ randn(elt, 3, 3)
   d = eig_vals(a)
   @test d ≈ diag(eig_full(a)[1])
 
-  d, v = eigh_full(ah)
-  @test ah * v ≈ v * d
+  a = Hermitian(randn(elt, 2, 2)) ⊗ Hermitian(randn(elt, 3, 3))
+  d, v = eigh_full(a)
+  @test a * v ≈ v * d
+
+  a = Hermitian(randn(elt, 2, 2)) ⊗ Hermitian(randn(elt, 3, 3))
+  @test_throws MethodError eigh_trunc(a)
+
+  a = Hermitian(randn(elt, 2, 2)) ⊗ Hermitian(randn(elt, 3, 3))
+  d = eigh_vals(a)
+  @test d ≈ diag(eigh_full(a)[1])
 
-  @test_throws MethodError eigh_trunc(ah)
+  a = randn(elt, 3, 2) ⊗ randn(elt, 4, 3)
+  n = left_null(a)
+  @test norm(n' * a) ≈ 0 atol = √eps(real(elt))
 
-  d = eigh_vals(ah)
-  @test d ≈ diag(eigh_full(ah)[1])
+  a = randn(elt, 2, 3) ⊗ randn(elt, 3, 4)
+  n = right_null(a)
+  @test norm(a * n') ≈ 0 atol = √eps(real(elt))
 end