Solving for Möbius transformation coefficients

Möbius transformations are functions of the form

$f(z) = \frac{az + b}{cz + d}$

where ad – bc ≠ 0.

A Möbius transformation is uniquely determined by its values at three points. Last year I wrote a post that mentioned how to determine the coefficients of a Möbius transformation. There I said

The unique bilinear transform sending z₁, z₂, and z₃ to w₁, w₂, and w₃ is given by

$\frac{(w - w_2)(w_3 - w_1)}{(w - w_1)(w_3 - w_2)} = \frac{(z - z_2)(z_3 - z_1)}{(z - z_1)(z_3 - z_2)}$

Plug in your constants and solve for w = f(z).

This is correct, but it still leaves a bit of work to do. In a particular case it’s not hard to find the coefficients, but it would be harder to find the coefficients in general.

There is an explicit formula for each of the parameters a, b, c, and d given the specified points z₁, z₂, z₃ and their images w₁, w₂, w₃ which we will present shortly. We could easily code-up these formulas except for one complication: we may want one of our inputs our outputs to be ∞. This is not merely an edge case: in applications, say to signal processing, you often want to specify the location of poles.

Simplest case: no infinities

If none of the inputs or outputs are infinite, the coefficients are given by

$\begin{align*} a &= \begin{vmatrix} z_1w_1 & w_1 & 1 \\ z_2w_2 & w_2 & 1 \\ z_3w_3 & w_3 & 1 \end{vmatrix} \\ b &= \begin{vmatrix} z_1w_1 & z_1 & w_1 \\ z_2w_2 & z_2 & w_2 \\ z_3w_3 & z_3 & w_3 \end{vmatrix} \\ c &= \begin{vmatrix} z_1\phantom{w_1} & w_1 & 1 \\ z_2\phantom{w_1} & w_2 & 1 \\ z_3\phantom{w_1} & w_3 & 1 \end{vmatrix} \\ d &= \begin{vmatrix} z_1w_1 & z_1 & 1 \\ z_2w_2 & z_2 & 1 \\ z_3w_3 & z_3 & 1 \end{vmatrix} \end{align*}$

This is one solution; multiplying all four coefficients by a non-zero constant gives another solution. But modulo a constant the solution is unique. To put it another way, the Möbius transformation is unique but it’s representation is only unique up to a constant multiplying the numerator and denominator.

You might hope for a second that this formula would work if you just let floating point infinities take care of themselves. But if any of the z‘s or w‘s is infinite, every determinant above is infinite.

Dealing with infinities

The possible cases of infinities are:

No infinities
Only one z is infinite
Only one w is infinite
A z and a w with the same subscript are infinite
A z and a w with different subscripts are infinite

The order of our z‘s and w‘s is arbitrary, so we can rearrange them if necessary for our convenience. So without loss of generality we may assume

No infinities
Only z₁ is infinite
Only w₁ is infinite
z₁ = w₁ = ∞
z₁ = w₂ = ∞

To handle the case (2) of z₁= ∞, divide each of the equations above by z₁and take the limit as z₁ approaches ∞. Since dividing one row of a matrix by a constant divides its determinant by the same amount, in each case we divide the first row by z₁. This works out nicely because z₁ only ever appears in the first row of each matrix. We can handle the case (3) of 2₁= ∞ analogously.

To handle the case (4) of z₁ = w₁ = ∞ we divide the first rows by z₁ w₁ and take the limit. To handle the case (5) of z₁ = w₂ = ∞ we divide the first rows by z₁ and the second rows by w₁ and take a limit.

If you find this limiting business dubious, it doesn’t matter: if the result is correct, it’s correct. And since Möbius transformations are determined by their values at three points, you can verify that each case is correct by sticking in z₁, z₂, z₃ and checking that you get w₁, w₂, w₃ out. You could do this manually, which I have, or trust the output of the code at the bottom of the post.

So here are our solutions.

(1) No infinities

Given above.

(2) Only z₁ = ∞

a = w₁ (w₂ – w₃)
b = w₁ (z₂ w₃ – z₃ w₂) + w₂ w₃ (z₃ – z₂)
c = w₂ – w₃
d = w₁ (z₂ – z₃) – z₂w₂ + z₃ w₃

(3) Only w₁ = ∞

a = z₁ (w₂ – w₃) – z₂w₂ + z₃ w₃
b = z₁ (z₂ w₃ – z₃ w₂) + z₂ z₃ (w₂ – w₃)
c = z₃ – z₂
d = z₁ (z₂ – z₃)

NB: You could derive these from the case z₁ = ∞ by inverting the transform. This means swapping z‘s with w‘s, swapping a and d, and negating b and c.

(4) z₁ = w₁ = ∞

a = w₂ – w₃
b = z₂ w₃ – z₃ w₂
c = 0
d = z₂ – z₃

(5) z₁ = w₂ = ∞

a = w₁
b = –z₂ w₃ + z₃ (w₃ – w₁)
c = 1
d = –z₂

Python code

import numpy as np

def all_finite(z1, z2, z3, w1, w2, w3):
    a = np.linalg.det(
            [[z1*w1, w1, 1],
             [z2*w2, w2, 1],
             [z3*w3, w3, 1]])
    b = np.linalg.det(
            [[z1*w1, z1, w1],
             [z2*w2, z2, w2],
             [z3*w3, z3, w3]])
    c = np.linalg.det(
            [[z1, w1, 1],
             [z2, w2, 1],
             [z3, w3, 1]])
    d = np.linalg.det(
            [[z1*w1, z1, 1],
             [z2*w2, z2, 1],
             [z3*w3, z3, 1]])
    return (a, b, c, d)

def z1_infinite(z1, z2, z3, w1, w2, w3):
    assert(np.isinf(z1))
    a = w1*(w2 - w3)
    b = w1*(z2*w3 - z3*w2) + w2*w3*(z3 - z2)
    c = w2 - w3
    d = w1*(z2 - z3) - z2*w2 + z3*w3
    return (a, b, c, d)    

def w1_infinite(z1, z2, z3, w1, w2, w3):
    assert(np.isinf(w1))
    a = z1*(w2 - w3) - z2*w2 + z3*w3
    b = z1*(z2*w3 - z3*w2) + z2*z3*(w2 - w3)
    c = z3 - z2
    d = z1*(z2 - z3)
    return (a, b, c, d)        

def z1w1_infinite(z1, z2, z3, w1, w2, w3):
    assert(np.isinf(z1) and np.isinf(w1))
    a = w2 - w3
    b = z2*w3 - z3*w2
    c = 0
    d = z2 - z3
    return (a, b, c, d)            

def z1w2_infinite(z1, z2, z3, w1, w2, w3):
    assert(np.isinf(z1) and np.isinf(w2))
    a = w1
    b = -z2*w3 + z3*(w3 - w1)
    c = 1
    d = -z2
    return (a, b, c, d)                

def mobius_coeff(z1, z2, z3, w1, w2, w3):

    infz = np.isinf(z1) or np.isinf(z2) or np.isinf(z3)
    infw = np.isinf(w1) or np.isinf(w2) or np.isinf(w3)    

    if infz:
        if np.isinf(z2):
            z1, z2 = z2, z1
            w1, w2 = w2, w1
        if np.isinf(z3):
            z1, z3 = z3, z1
            w1, w3 = w3, w1
        if infw:
            if np.isinf(w1):
                return z1w1_infinite(z1, z2, z3, w1, w2, w3)
            if np.isinf(w3):
                z2, z3 = z3, z2
                w2, w3 = w3, w2
            return z1w2_infinite(z1, z2, z3, w1, w2, w3)
        else:
            return z1_infinite(z1, z2, z3, w1, w2, w3)
        
    if infw: # and all z finite
        if np.isinf(w2):
            z1, z2 = z2, z1
            w1, w2 = w2, w1
        if np.isinf(w3):
            z1, z3 = z3, z1
            w1, w3 = w3, w1
        return w1_infinite(z1, z2, z3, w1, w2, w3)

    return all_finite(z1, z2, z3, w1, w2, w3)

def mobius(x, a, b, c, d):
    if np.isinf(x):
        if c == 0:
            return np.inf
        return a/c
    if c*x + d == 0:
        return np.inf
    else:
        return (a*x + b)/(c*x + d)
    
def test_mobius(z1, z2, z3, w1, w2, w3):
    tolerance = 1e-6
    a, b, c, d = mobius_coeff(z1, z2, z3, w1, w2, w3)
    for (x, y) in [(z1, w1), (z2, w2), (z3, w3)]:
        m = mobius(x, a, b, c, d)
        assert(np.isinf(m) and np.isinf(y) or abs(m - y) <= tolerance)

test_mobius(1, 2, 3, 6, 4, 2)
test_mobius(1, 2j, 3+7j, 6j, -4, 2) 
test_mobius(np.inf, 2, 3, 8j, -2, 0)
test_mobius(0, np.inf, 2, 3, 8j, -2)
test_mobius(0, -1, np.inf, 2, 8j, -2)
test_mobius(1, 2, 3, np.inf, 44j, 0)
test_mobius(1, 2, 3, 1, np.inf, 40j)
test_mobius(-1, 0, 3j, 1, -1j, np.inf)
test_mobius(np.inf, -1j, 5, np.inf, 2, 8)
test_mobius(1, np.inf, -1j, 5, np.inf, 2)
test_mobius(12, 0, np.inf, -1j, 5, np.inf)
test_mobius(np.inf, -1j, 5, 0, np.inf, -1)
test_mobius(6, np.inf, -1j, 0, 8, np.inf)
test_mobius(6, 3j, np.inf, -1j, np.inf, 1)

More Möbius transformation posts

2 thoughts on “Solving for Möbius transformation coefficients”

Jan Van lent

28 March 2021 at 20:19

Using the scipy.linalg.null_space command can make this a slightly more concise.
Each pair of points corresponds to an equation of the form
z a + b – z w c – w d = 0.
We divide through by z and/or w if they are infinite.
I did have to relax the test, so it does not give exactly the same kind of solution.
Not sure if posting code works.

import numpy as np
import scipy.linalg as sla

# w = (a*z+b)/(c*z+d)
# z a + b - z w c - w d = 0

def mobius_row(z, w):
    if np.isinf(z):
        if np.isinf(w):
            row = [0, 0, -1, 0]
        else:
            row = [1, 0, -w, 0]
    else:
        if np.isinf(w):
            row = [0, 0, -z, -1]
        else:
            row = [z, 1, -z*w, -w]
    return row

def mobius_coeff(z1, z2, z3, w1, w2, w3):
    A = np.array([ mobius_row(z1, w1),
                   mobius_row(z2, w2),
                   mobius_row(z3, w3) ])
    a, b, c, d = sla.null_space(A).flat
    return a, b, c, d

assert(np.isinf(m) and np.isinf(y) or
       np.isinf(m) and 1/abs(y) <= tolerance or
       np.isinf(y) and 1/abs(m) <= tolerance or
       abs(m - y) <= tolerance)

John

28 March 2021 at 20:29

Thanks. That’s neat.

I edited your comment to turn <code> into <pre> and that preserved your code formatting.

Comments are closed.