Finding linear least square curve fit with restrictions

Let's assume that we have a set of eight different lights each with their own spectrum (image 1). Each light can be separattely dimmed from 0 to 100% as shown in the image 2.

Image 1

Image 2

Image 3

Image 3 shows four example conditions that we want to match with our lighting setup. The blue dashed line is the target and the green solid line is the best match.

The best match is acquired using dimming levels:

Table 1

How do we get there values?


Mathematical setup

Let's model the problem with matrices.

The matrix A is the information we have about the lights. The matrix A remains always constant. The size of A is 85 rows and 8 columns. Each column represents a single light and each row provides information about particular light frequency.

The vector t represents the target condition that we want to achieve. The size of t is 85 rows and 1 column.

The vector x is the unknown which we want to calculate.

The vector b is the light condition that we actually achieve.

Euclidean norm provides a base to compare different estimations. The euclidean norm takes the square root of the sum of individual errors to the second power.

Moore-Penrose pseudoinverse is proven to provide the solution if certain conditions are met - in this case the lights need to be different from each other in a way that no single light can be simulated using the other lights.


Solution

Excel 2013

Excel file

MMULT "Note: To work correctly, the formula in the example needs to be entered as an array formula in the Excel program. After copying the example to a blank worksheet, select the range C8:D9 starting with the formula cell. Press F2, and then press Ctrl+Shift+Enter. If the formula is not entered as an array formula, a single result (2) will be returned in cell C8."

// Matrix A
A2:D86

// transpose(A)
{=TRANSPOSE(A2:D86)}

// transpose(A) * A
{=MMULT(A89:CG92;A2:D86)}

// inverse(transpose(A) * A)
{=MINVERSE(A95:D98)}

// pseudoinverse(A) = inverse(transpose(A) * A) * transpose(A)
{=MMULT(A101:D104;A89:CG92)}
Problem with pseudoinverse

Image 3 has a different target (blue line). Now the best fit using pseudoinverse (red line) shows a problem. The optimal dimming levels for light 1 and 4 are negative. The dimming levels cannot be below 0% or more than 100%.

Image 4

Image 5

2b (yellow line) uses bounded values:

=MAX(MIN(A2;1);0)
=MAX(MIN(A3;1);0)
=MAX(MIN(A4;1);0)
=MAX(MIN(A5;1);0)

2c (green line) uses Excel GRG-nonlinear solver

The solver target is to minimize the error.

Node.js

Node.js with math.js provides a base for the implementation.

pseudoInverse.js
/**************************************************************/
/*
 * pseudoinverse.js
 *                                        *                      
 * Calculate Moore-Penrose pseudoinverse A
 * http://en.m.wikipedia.org/wiki/Moore-Penrose_pseudoinverse
 *
 *    *     T  -1  T
 *   A  = (A A)   A
 *
 */

var math = require('mathjs');

module.exports = {
    pseudoInverse: function (matrix) {

        var type = matrix.type;
        if (type != 'DenseMatrix') {
            throw new Exception('Invalid type');
        }
        var size = matrix.size();
        if (size.length != 2) {
            throw new Exception('Invalid size');
        }
        if (size[1] > size[0]) {
            throw new Exception('Invalid size');
        }

        var matrixTranspose = math.transpose(matrix);
        var matrixMultiply = math.multiply(matrixTranspose,
            matrix);
        if (math.det(matrixMultiply) == 0) {
            throw new Exception('Invalid matrix');
        }
        return math.multiply(math.inv(matrixMultiply),
                            matrixTranspose);   
    }
}   
unboundedLinearLeastSquaresCurveFit.js
/**************************************************************/
/* unboundedLinearLeastSquaresCurveFit.js                         
 * Calculate best fit x for Ax = b    * 
 * using Moore-Penrose pseudoinverse A
 * 
 *
 *   Ax = b
 *
 *        *
 *   x = A b
 *
 *    *     T  -1  T
 *   A  = (A A)   A
 *
 */

var math = require('mathjs');
var pInv = require('./pseudoInverse');
var mxA = require('./createMatrix'); 

module.exports = {
    unboundedLinearLeastSquaresCurveFit: 
        function (vectors, target, forceMatrix) {
            var A = mxA.createMatrix(vectors);
        return this.unboundedLinearLeastSquaresCurveFitA(A,
            target);
    },
    unboundedLinearLeastSquaresCurveFitA: 
        function (A, target, forceMatrix) {
            var x = math.multiply(pInv.pseudoInverse(A), target);
            if (forceMatrix && typeof (x) == 'number') { 
                // convert number back to 1x1 matrix
                return math.matrix([[x]]);
        }
        else return x;
    }
}

Notes

N-vector is different than N x 1-matrix

...
var b = math.matrix(target);
b.resize([target.length, 1]);
...

1x1-matrix is different than number (scalar)

...
var x = math.multiply(pA, b);
if (typeof (x) == 'number') { // convert number back to 1x1 matrix
    return math.matrix([[x]]);
}
else return x;
...

Github repository

NPM package