finding matrix through optimisation

I am looking for algorithm to solve the following problem :

I have two sets of vectors, and I want to find the matrix that best approximate the transformation from the inpu开发者_如何学Pythont vectors to the output vectors.

vectors are 3x1, so matrix is 3x3.

This is the general problem. My particular problem is I have a set of RGB colors, and another set that contains the desired color. I am trying to find an RGB to RGB transformation that would give me colors closer to the desired ones.

There is correspondence between the input and output vectors, so computing an error function that should be minimized is the easy part. But how can I minimize this function ?

---

This is a classic linear algebra problem, the key phrase to search on is "multiple linear regression".

I've had to code some variation of this many times over the years. For example, code to calibrate a digitizer tablet or stylus touch-screen uses the same math.

---

Here's the math:

Let p be an input vector and q the corresponding output vector.

The transformation you want is a 3x3 matrix; call it A.

For a single input and output vector p and q, there is an error vector e

e = q - A x p

The square of the magnitude of the error is a scalar value:

eT x e = (q - A x p)T x (q - A x p)

(where the T operator is transpose).

What you really want to minimize is the sum of e values over the sets:

E = sum (e)

This minimum satisfies the matrix equation D = 0 where

D(i,j) = the partial derivative of E with respect to A(i,j)

Say you have N input and output vectors.

Your set of input 3-vectors is a 3xN matrix; call this matrix P. The ith column of P is the ith input vector.

So is the set of output 3-vectors; call this matrix Q.

When you grind thru all of the algebra, the solution is

A = Q x PT x (P x PT) ^-1

(where ^-1 is the inverse operator -- sorry about no superscripts or subscripts)

---

Here's the algorithm:

Create the 3xN matrix P from the set of input vectors.

Create the 3xN matrix Q from the set of output vectors.

Matrix Multiply R = P x transpose (P)

Compute the inverseof R

Matrix Multiply A = Q x transpose(P) x inverse (R)

using the matrix multiplication and matrix inversion routines of your linear algebra library of choice.

---

However, a 3x3 affine transform matrix is capable of scaling and rotating the input vectors, but not doing any translation! This might not be general enough for your problem. It's usually a good idea to append a "1" on the end of each of the 3-vectors to make then a 4-vector, and look for the best 3x4 transform matrix that minimizes the error. This can't hurt; it can only lead to a better fit of the data.

---

You don't specify a language, but here's how I would approach the problem in Matlab.

• v1 is a 3xn matrix, containing your input colors in vertical vectors
• v2 is also a 3xn matrix containing your output colors

You want to solve the system

M*v1 = v2
M = v2*inv(v1)

However, v1 is not directly invertible, since it's not a square matrix. Matlab will solve this automatically with the mrdivide operation (M = v2/v1), where M is the best fit solution.

eg: 
>> v1 = rand(3,10);
>> M = rand(3,3);
>> v2 = M * v1;
>> v2/v1 - M

ans =

   1.0e-15 *

    0.4510    0.4441   -0.5551
    0.2220    0.1388   -0.3331
    0.4441    0.2220   -0.4441

>> (v2 + randn(size(v2))*0.1)/v1 - M
ans =

    0.0598   -0.1961    0.0931
   -0.1684    0.0509    0.1465
   -0.0931   -0.0009    0.0213

This gives a more language-agnostic solution on how to solve the problem.

---

Some linear algebra should be enough :

Write the average squared difference between inputs and outputs ( the sum of the squares of each difference between each input and output value ). I assume this as definition of "best approximate"

This is a quadratic function of your 9 unknown matrix coefficients.

To minimize it, derive it with respect to each of them.

You will get a linear system of 9 equations you have to solve to get the solution ( unique or a space variety depending on the input set )

When the difference function is not quadratic, you can do the same but you have to use an iterative method to solve the equation system.