projected move regression

gistfile1.jl

using LinearAlgebra

function linreg(xs, ys)
    # Axs + b ≈ ys
    @assert size(xs,2) == size(ys,2)
    nobs = size(xs, 2)
    xs_ = [xs; transpose(ones(nobs))]
    ys_ = [ys; transpose(ones(nobs))]
    A_ = ys_ / xs_
    @assert A_[end, end] ≈ 1
    @assert norm(A_[end, 1:end-1]) < sqrt(eps(eltype(A_)))
    A = A_[1:end-1, 1:end-1]
    b = A_[1:end-1, end]
    A, b
end

function projected_move_regression(proj, xs, ys)
    # ys ≈ proj * R * xs .+ b
    # where R is a rotation
    projR, b = linreg(xs, ys)
    
    A = proj \ projR
    U,S,V = svd(A)
    Vt = transpose(V)
    @assert A ≈ U * Diagonal(S) * Vt
    R = U * Vt
    R, b
end


using Rotations, StaticArrays

using Test
A = randn(2,2)
b = randn(2)
xs = randn(2,10)
ys = A * xs .+ b 

A_reg, b_reg = linreg(xs, ys)
@test A_reg ≈ A
@test b_reg ≈ b

rot_truth = rand(RotMatrix{3}) 
proj_trans_truth = @SVector randn(2)
proj = [1 0 0; 0 1 0]
N = 5
xs = randn(3,N)
ys = proj * rot_truth * xs .+ proj_trans_truth

R, b = projected_move_regression(proj, xs, ys)
@test R ≈ rot_truth
@test b ≈ proj_trans_truth