Como enviar cartas pela internet

November 27, 2009 / cesarsouza / 1 Comment

Pode parecer estranho neste mundo onde reina a internet e o correio digital. Mas ainda assim há horas em que precisamos recorrer ao ~~bom e~~ velho correio convencional para enviar correspondência. Além do que, uma carta em papel é muito mais formal e impactante do que um email. No entanto, se você está com preguiça de ir até o correio apenas para postar uma carta, ou se precisa enviar correspondência para outro país e não está certo de como preencher o envelope nos padrões do país destinatário, não se preocupe! Existem serviços que permitem enviar cartas convencionais através da internet de maneira eficiente e confiável!

Para enviar cartas pela internet, existem vários sites que podemos utilizar.

Para o território nacional

Serviço próprio dos correios. Uma carta custa, em média, R$3,60. Pode ser usada para cartas internacionais também, mas é um tanto caro e as opções de pagamento incluem penas cartão de crédito.

http://www.correios.com.br/produtos_servicos/catalogo/carta_internet.cfm

Para outros países

O melhor serviço que encontrei (e já utilizei) é o Mail A Letter. Uma carta internacional custa U$1,99 (R$3,53 em out/2009) e sua maior vantagem é que o site aceita pagamentos via paypal. Para cartas para os Estados Unidos, o serviço custa apenas U$0,99.

http://www.mailaletter.com/

Neural Network Learning by the Levenberg-Marquardt Algorithm with Bayesian Regularization (part 1)

November 18, 2009 / cesarsouza / 58 Comments

A complete explanation for the totally lost, part 1 of 2.

Downloads

Download the sample applications.
Download the Accord.NET Framework.
Browse the source code online.

The code has been incorporated in Accord.NET Framework, which includes the latest version of this code plus many other statistics and machine learning tools.

Overview

The problem of neural network learning can be seen as a function optimization problem, where we are trying to determine the best network parameters (weights and biases) in order to minimize network error. This said, several function optimization techniques from numerical linear algebra can be directly applied to network learning, one of these techniques being the Levenberg-Marquardt algorithm.
The Levenberg–Marquardt algorithm provides a numerical solution to the problem of minimizing a (generally nonlinear) function, over a space of parameters for the function. It is a popular alternative to the Gauss-Newton method of finding the minimum of a function.

Neural Networks

Neural networks are a relatively new artificial intelligence technique. In most cases an ANN is an adaptive system that changes its structure based on external or internal information that flows through the network during the learning phase. The learning procedure tries is to find a set of connections w that gives a mapping that fits the training set well.
Furthermore, neural networks can be viewed as highly nonlinear functions with the basic the form:

F(x,w) = y

(1)

Where x is the input vector presented to the network, w are the weights of the network, and y is the corresponding output vector approximated or predicted by the network. The weight vector w is commonly ordered first by layer, then by neurons, and finally by the weights of each neuron plus its bias.
This view of network as an parameterized function will be the basis for applying standard function optimization methods to solve the problem of neural network training.

AForge Framework

AForge.NET Framework is a C# framework designed for developers and researchers in the fields of Computer Vision and Artificial Intelligence. Here, the Levenberg-Marquardt learning algorithm is implemented as a class implementing the ISupervisedLearning interface from the AForge framework.

Network training as a function optimization problem

As mentioned previously, neural networks can be viewed as highly non-linear functions. From this perspective, the training problem can be considered as a general function optimization problem, with the adjustable parameters being the weights and biases of the network, and the Levenberg-Marquardt can be straightforward applied in this case.

Levenberg-Marquardt Algorithm

The Levenberg-Marquardt algorithm is a very simple, but robust, method for approximating a function. Basically, it consists in solving the equation:

(J^tJ + λI)δ = J^tE

(2)

Where J is the Jacobian matrix for the system, λ is the Levenberg’s damping factor, δ is the weight update vector that we want to find and E is the error vector containing the output errors for each input vector used on training the network. The δ tell us by how much we should change our network weights to achieve a (possibly) better solution. The J^tJ matrix can also be known as the approximated Hessian.
The λ damping factor is adjusted at each iteration, and guides the optimization process. If reduction of E is rapid, a smaller value can be used, bringing the algorithm closer to the Gauss–Newton algorithm, whereas if an iteration gives insufficient reduction in the residual, λ can be increased, giving a step closer to the gradient descent direction.

Computing the Jacobian

The Jacobian is a matrix of all first-order partial derivatives of a vector-valued function. In the neural network case, it is a N-by-W matrix, where N is the number of entries in our training set and W is the total number of parameters (weights + biases) of our network. It can be created by taking the partial derivatives of each output in respect to each weight, and has the form:

Where F(x_i, w) is the network function evaluated for the i-th input vector of the training set using the weight vector w and w_j is the j-th element of the weight vector w of the network.
In traditional Levenberg-Marquardt implementations, the Jacobian is approximated by using finite differences. However, for neural networks, it can be computed very efficiently by using the chain rule of calculus and the first derivatives of the activation functions.

Approximating the Hessian

For the least-squares problem, the Hessian generally doesn’t needs to be calculated. As stated earlier, it can be approximated by using the Jacobian matrix with the formula:

H ≈ J^tJ

(3)

Which is is a very good approximation of the Hessian if the residual errors at the solution are “small”. If the residuals are not sufficiently small at the solution, this approach may result in slow convergence. The Hessian can also be used to apply regularization to the learning process, which will be discussed later.

Solving the Levenberg-Marquardt equation

Levenberg’s main contribution to the method was the introduction of the damping factor λ. This value is summed to every member of the approximate Hessian diagonal before the system is solved for the gradient. Tipically, λ would start as a small value such as 0.1.
Then, the Levenberg-Marquardt equation is solved, commonly by using a LU decomposition. However, the system can only be solved if the approximated Hessian has not become singular (not having an inverse). If this is the case, the equation can still be solved by using a SVD decomposition.
After the equation is solved, the weights w are updated using δ and network errors for each entry in the training set are recalculated. If the new sum of squared errors has decreased, λ is decreased and the iteration ends. If it has not, then the new weights are discarded and the method is repeated with a higher value for λ.
This adjustment for λ is done by using an adjustment factor v, usually defined as 10. If λ needs to increase, it is multiplied by v. If it needs to decrease, then it is divided by v. The process is repeated until the error decreases. When this happens, the current iteration ends.

General Levenberg-Marquardt Algorithm

As stated earlier, the Levenberg-Marquardt consists basically in solving (2) with different λ values until the sum of squared error decreases. So, each learning iteration (epoch) will consist of the following basic steps:

Compute the Jacobian (by using finite differences or the chain rule)

Compute the error gradient

g = J^tE

Approximate the Hessian using the cross product Jacobian (eq. 3)

H = J^tJ

Solve (H + λI)δ = g to find δ

Update the network weights w using δ

Recalculate the sum of squared errors using the updated weights

If the sum of squared errors has not decreased,

Discard the new weights, increase λ using v and go to step 4.

Else decrease λ using v and stop.

Variations of the algorithm may include different values for v, one for decreasing λ and other for increasing it. Others may solve (H + λdiag(H))δ = g instead of (H + λI)δ = g (2), while others may select the initial λ according to the size of the elements on H, by setting λ₀ = t max(diag(H)), where t is a value chosen by the user. I’ve chosen the identity matrix equation because, apparently, it is the same method implemented internally by the Neural Network Toolbox in MATLAB.
We can see we will have a problem if the error does not decrease after a some iterations. In this case, the algorithm also stops if λ becomes too large.

Limitations

The Levenberg-Marquardt is very sensitive to the initial network weighs. Also, it does not consider outliers in the data, what may lead to overfitting noise. To avoid those situations, we can use a technique known as regularization.

In the next part of this article (part 2), we’ll discuss more about Bayesian regularization. We will also present and demonstrate the usage of the article’s accompanying source code. Please click here to go to Neural Network Learning by the Levenberg-Marquardt Algorithm with Bayesian Regularization (part 2).

Go to Neural Network Learning by the Levenberg-Marquardt Algorithm (part 2)

CodeProject

Neural Network Learning by the Levenberg-Marquardt Algorithm with Bayesian Regularization (part 2)

November 18, 2009 / cesarsouza / 49 Comments

A complete explanation for the totally lost, part 2 of 2. Click here to return to part 1.

Downloads

Download the sample applications.
Download the Accord.NET Framework.
Browse the source code online.

The code has been incorporated in Accord.NET Framework, which includes the latest version of this code plus many other statistics and machine learning tools.

Bayesian Regularization

To overcome the problem in interpolating noisy data, MacKay (1992) has proposed a Bayesian framework which can be directly applied to the neural network learning problem. It also allows to estimate the effective number of parameters actually used by the model – in this case, the number of network weights actually needed to solve a particular problem.
Bayesian regularization expands the cost function to search not only for the minimal error, but for the minimal error using the minimal weights. It works by introducing two Bayesian hyperparameters, alpha and beta, to tell which direction (minimal error or minimal weights) the learning process must seek.
The cost function will then become:

C(k) = β*E_d+α*E_w, where:

E_d is the sum of squared errors, and

E_w is the sum of squared weights

By using Bayesian regularization, one can avoid costly cross validation. It is particularly useful to problems that can not, or would suffer, if a portion of the available data were reserved to a validation set. Regularization also reduces (or eliminates) the need for testing different number of hidden neurons for a problem. A third variable, gamma, indicates the number of effective weights being used by the network, thus giving an indication on how complex the network should be.
Many practical realizations of Bayesian Regularization perform an update of the hyperparameters alpha and beta after each training cycle. However, according to Poland (2001) the most popular update algorithm fails to produce robust iterates if there is not much training data. The following algorithm shows how to update the hyperparameters according to both rules. Both methods base their calculations in the inverse Hessian matrix.

Expanded Levenberg-Marquardt Algorithm

Adding bayesian regularization to the Levenberg-Marquardt adds little overhead to the process because we already have an Hessian approximation. So, the algorithm for a learning epoch then becomes:

Compute the Jacobian (by finite differences or using the chain rule)

Compute the error gradient

g = J^tE

Approximate the Hessian using the cross product Jacobian

H = J^tJ

Calculate the cost function

C = β*E_d+α*E_w, where:

E_d is the sum of squared errors, and

E_w is the sum of squared weights

Solve (H + λI)δ = g to find δ

Update the network weights w using δ

Recalculate the cost function using the updated weights

If the cost has not decreased,

Discard the new weights, increase λ using v and go to step 5.

Else decrease λ using v

Update the Bayesian hyperparameters using MacKay’s or Poland’s formulae:

gamma = W – (alpha * tr(H^-1))

beta = (N – gamma) / 2.0 * E_d

alpha = W / (2.0 * E_w + tr(H-1)) [modified Poland’s update], or

alpha = gamma / (2.0 * E_w) [original MacKay’s update], where:

W is the number of network parameters (number of weights and biases)

N is the number of entries in the training set

tr(H^-1) is the trace of the inverse Hessian matrix

Source Code

Here is the code for running a learning epoch of the Levenberg-Marquardt algorithm.

<span>/// &lt;summary&gt;</span>
<span>///   Runs learning epoch.</span>
<span>/// &lt;/summary&gt;</span>
<span>/// </span>
<span>/// &lt;param name="input"&gt;Array of input vectors.&lt;/param&gt;</span>
<span>/// &lt;param name="output"&gt;Array of output vectors.&lt;/param&gt;</span>
<span>/// </span>
<span>/// &lt;returns&gt;Returns summary learning error for the epoch.&lt;/returns&gt;</span>
<span>/// </span>
<span>/// &lt;remarks&gt;&lt;para&gt;The method runs one learning epoch, by calling running necessary</span>
<span>/// iterations of the Levenberg Marquardt to achieve an error decrease.&lt;/remarks&gt;</span>
<span>///</span>
<span>public</span> <span>double</span> RunEpoch(<span>double</span>[][] input, <span>double</span>[][] output)
{
    <span>// Initial definitions and memmory allocations</span>
    <span>int</span> N = input.Length;
    <span>this</span>.errors = <span>new</span> <span>double</span>[N];
    <span>this</span>.jacobian = <span>new</span> <span>double</span>[N, numberOfParameters];
    LuDecomposition decomposition = <span>null</span>;
    <span>double</span> sumOfSquaredErrors = 0.0;
    <span>double</span> sumOfSquaredWeights = 0.0;
    <span>double</span> trace = 0.0;


    <span>// Compute the Jacobian matrix</span>
    <span>if</span> (method == JacobianMethod.ByBackpropagation)
    {
        sumOfSquaredErrors = JacobianByChainRule(input, output);
    }
    <span>else</span>
    {
        sumOfSquaredErrors = JacobianByFiniteDifference(input, output);
    }


    <span>// Create the initial weights vector</span>
    sumOfSquaredWeights = CreateWeights();


    <span>// Compute Jacobian errors and Hessian</span>
    <span>for</span> (<span>int</span> i = 0; i &lt; numberOfParameters; i++)
    {
        <span>// Compute Jacobian Matrix Errors</span>
        <span>double</span> s = 0.0;
        <span>for</span> (<span>int</span> j = 0; j &lt; N; j++)
        {
            s += jacobian[j, i] * errors[j];
        }
        gradient[i] = s;

        <span>// Compute Quasi-Hessian Matrix using Jacobian (H = J'J)</span>
        <span>for</span> (<span>int</span> j = 0; j &lt; numberOfParameters; j++)
        {
            <span>double</span> c = 0.0;
            <span>for</span> (<span>int</span> k = 0; k &lt; N; k++)
            {
                c += jacobian[k, i] * jacobian[k, j];
            }
            hessian[i, j] = beta * c;
        }
    }

    <span>// Store the Hessian diagonal for future computations</span>
    <span>for</span> (<span>int</span> i = 0; i &lt; numberOfParameters; i++)
        diagonal[i] = hessian[i, i];


    <span>// Define the objective function</span>
    <span>// bayesian regularization objective function</span>
    <span>double</span> objective = beta * sumOfSquaredErrors + alpha * sumOfSquaredWeights;
    <span>double</span> current = objective + 1.0;



    <span>// Begin of the main Levenberg-Macquardt method</span>
    lambda /= v;

    <span>// We'll try to find a direction with less error</span>
    <span>//  (or where the objective function is smaller)</span>
    <span>while</span> (current &gt;= objective &amp;&amp; lambda &lt; lambdaMax)
    {
        lambda *= v;

        <span>// Update diagonal (Levenberg-Marquardt formula)</span>
        <span>for</span> (<span>int</span> i = 0; i &lt; numberOfParameters; i++)
            hessian[i, i] = diagonal[i] + (lambda + alpha);

        <span>// Decompose to solve the linear system</span>
        decomposition = <span>new</span> LuDecomposition(hessian);

        <span>// Check if the Jacobian has become non-invertible</span>
        <span>if</span> (!decomposition.NonSingular) <span>continue</span>;

        <span>// Solve using LU (or SVD) decomposition</span>
        deltas = decomposition.Solve(gradient);

        <span>// Update weights using the calculated deltas</span>
        sumOfSquaredWeights = UpdateWeights();

        <span>// Calculate the new error</span>
        sumOfSquaredErrors = 0.0;
        <span>for</span> (<span>int</span> i = 0; i &lt; N; i++)
        {
            network.Compute(input[i]);
            sumOfSquaredErrors += CalculateError(output[i]);
        }

        <span>// Update the objective function</span>
        current = beta * sumOfSquaredErrors + alpha * sumOfSquaredWeights;

        <span>// If the object function is bigger than before, the method</span>
        <span>//  is tried again using a greater dumping factor.</span>
    }

    <span>// If this iteration caused a error drop, then next iteration</span>
    <span>//  will use a smaller damping factor.</span>
    lambda /= v;



    <span>// If we are using bayesian regularization, we need to</span>
    <span>//   update the bayesian hyperparameters alpha and beta</span>
    <span>if</span> (useBayesianRegularization)
    {
        <span>// References: </span>
        <span>// - http://www-alg.ist.hokudai.ac.jp/~jan/alpha.pdf</span>
        <span>// - http://www.inference.phy.cam.ac.uk/mackay/Bayes_FAQ.html</span>

        <span>// Compute the trace for the inverse hessian</span>
        trace = Matrix.Trace(decomposition.Inverse());

        <span>// Poland update's formula:</span>
        gamma = numberOfParameters - (alpha * trace);
        alpha = numberOfParameters / (2.0 * sumOfSquaredWeights + trace);
        beta = System.Math.Abs((N - gamma) / (2.0 * sumOfSquaredErrors));
        <span>//beta = (N - gama) / (2.0 * sumOfSquaredErrors);</span>

        <span>// Original MacKay's update formula:</span>
        <span>//  gama = (double)networkParameters - (alpha * trace);</span>
        <span>//  alpha = gama / (2.0 * sumOfSquaredWeights);</span>
        <span>//  beta = (gama - N) / (2.0 * sumOfSquaredErrors);</span>
    }

    <span>return</span> sumOfSquaredErrors;
}

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

/// <summary>

/// Runs learning epoch.

/// </summary>

///

/// <param name="input">Array of input vectors.</param>

/// <param name="output">Array of output vectors.</param>

///

/// <returns>Returns summary learning error for the epoch.</returns>

///

/// <remarks><para>The method runs one learning epoch, by calling running necessary

/// iterations of the Levenberg Marquardt to achieve an error decrease.</remarks>

///

public double RunEpoch(double[][] input, double[][] output)

{

// Initial definitions and memmory allocations

int N = input.Length;

this.errors = new double[N];

this.jacobian = new double[N, numberOfParameters];

LuDecomposition decomposition = null;

double sumOfSquaredErrors = 0.0;

double sumOfSquaredWeights = 0.0;

double trace = 0.0;

// Compute the Jacobian matrix

if (method == JacobianMethod.ByBackpropagation)

{

sumOfSquaredErrors = JacobianByChainRule(input, output);

}

else

{

sumOfSquaredErrors = JacobianByFiniteDifference(input, output);

}

// Create the initial weights vector

sumOfSquaredWeights = CreateWeights();

// Compute Jacobian errors and Hessian

for (int i = 0; i < numberOfParameters; i++)

{

// Compute Jacobian Matrix Errors

double s = 0.0;

for (int j = 0; j < N; j++)

{

s += jacobian[j, i] * errors[j];

}

gradient[i] = s;

// Compute Quasi-Hessian Matrix using Jacobian (H = J'J)

for (int j = 0; j < numberOfParameters; j++)

{

double c = 0.0;

for (int k = 0; k < N; k++)

{

c += jacobian[k, i] * jacobian[k, j];

}

hessian[i, j] = beta * c;

}

// Store the Hessian diagonal for future computations

for (int i = 0; i < numberOfParameters; i++)

diagonal[i] = hessian[i, i];

// Define the objective function

// bayesian regularization objective function

double objective = beta * sumOfSquaredErrors + alpha * sumOfSquaredWeights;

double current = objective + 1.0;

// Begin of the main Levenberg-Macquardt method

lambda /= v;

// We'll try to find a direction with less error

// (or where the objective function is smaller)

while (current >= objective && lambda < lambdaMax)

{

lambda *= v;

// Update diagonal (Levenberg-Marquardt formula)

for (int i = 0; i < numberOfParameters; i++)

hessian[i, i] = diagonal[i] + (lambda + alpha);

// Decompose to solve the linear system

decomposition = new LuDecomposition(hessian);

// Check if the Jacobian has become non-invertible

if (!decomposition.NonSingular) continue;

// Solve using LU (or SVD) decomposition

deltas = decomposition.Solve(gradient);

// Update weights using the calculated deltas

sumOfSquaredWeights = UpdateWeights();

// Calculate the new error

sumOfSquaredErrors = 0.0;

for (int i = 0; i < N; i++)

{

network.Compute(input[i]);

sumOfSquaredErrors += CalculateError(output[i]);

}

// Update the objective function

current = beta * sumOfSquaredErrors + alpha * sumOfSquaredWeights;

// If the object function is bigger than before, the method

// is tried again using a greater dumping factor.

}

// If this iteration caused a error drop, then next iteration

// will use a smaller damping factor.

lambda /= v;

// If we are using bayesian regularization, we need to

// update the bayesian hyperparameters alpha and beta

if (useBayesianRegularization)

{

// References:

// - http://www-alg.ist.hokudai.ac.jp/~jan/alpha.pdf

// - http://www.inference.phy.cam.ac.uk/mackay/Bayes_FAQ.html

// Compute the trace for the inverse hessian

trace = Matrix.Trace(decomposition.Inverse());

// Poland update's formula:

gamma = numberOfParameters - (alpha * trace);

alpha = numberOfParameters / (2.0 * sumOfSquaredWeights + trace);

beta = System.Math.Abs((N - gamma) / (2.0 * sumOfSquaredErrors));

//beta = (N - gama) / (2.0 * sumOfSquaredErrors);

// Original MacKay's update formula:

// gama = (double)networkParameters - (alpha * trace);

// alpha = gama / (2.0 * sumOfSquaredWeights);

// beta = (gama - N) / (2.0 * sumOfSquaredErrors);

}

return sumOfSquaredErrors;

}

Using the Code

Code usage is very simple. It follows the common AForge.NET learning algorithm usage, instantiating a LevenbergMarquardtLearning class instead of the more traditional BackpropagationLearning. Optionally, we may opt for using the Jacobian approximated by finite differences instead of computing it directly (although it would not be very practical unless for debugging purposes). We may also enable bayesian regularization by passing a boolean parameter to the constructor.

<span>// initialize input and output values</span>
<span>double</span>[][] input = <span>new</span> <span>double</span>[4][] {
    <span>new</span> <span>double</span>[] {0, 0}, <span>new</span> <span>double</span>[] {0, 1},
    <span>new</span> <span>double</span>[] {1, 0}, <span>new</span> <span>double</span>[] {1, 1}
};
<span>double</span>[][] output = <span>new</span> <span>double</span>[4][] {
    <span>new</span> <span>double</span>[] {0}, <span>new</span> <span>double</span>[] {1},
    <span>new</span> <span>double</span>[] {1}, <span>new</span> <span>double</span>[] {0}
};
<span>// create neural network</span>
ActivationNetwork   network = <span>new</span> ActivationNetwork(
    SigmoidFunction( 2 ),
    2, <span>// two inputs in the network</span>
    2, <span>// two neurons in the first layer</span>
    1 ); <span>// one neuron in the second layer</span>
<span>// create teacher</span>
LevenbergMarquardtLearning teacher = <span>new</span> LevenbergMarquardtLearning( network );
<span>// loop</span>
<span>while</span> ( !needToStop )
{
    <span>// run epoch of learning procedure</span>
    <span>double</span> error = teacher.RunEpoch( input, output );
    <span>// check error value to see if we need to stop</span>
    <span>// ...</span>
}

// initialize input and output values

double[][] input = new double[4][] {

new double[] {0, 0}, new double[] {0, 1},

new double[] {1, 0}, new double[] {1, 1}

};

double[][] output = new double[4][] {

new double[] {0}, new double[] {1},

new double[] {1}, new double[] {0}

};

// create neural network

ActivationNetwork network = new ActivationNetwork(

SigmoidFunction( 2 ),

2, // two inputs in the network

2, // two neurons in the first layer

1 ); // one neuron in the second layer

// create teacher

LevenbergMarquardtLearning teacher = new LevenbergMarquardtLearning( network );

// loop

while ( !needToStop )

{

// run epoch of learning procedure

double error = teacher.RunEpoch( input, output );

// check error value to see if we need to stop

// ...

}

Sample applications

The original sample applications from AForge.NET have been modified to use Levenberg-Marquardt learning. However, it is quite tricky to get them to converge adequately.

Function approximation using regular Levenberg-Marquardt

Function approximation using the Levenberg-Marquardt with Bayesian regularization. Note that some of them does not converge to a best fit solution; instead, a possibly more general solution is preferred. The top-right image uses the Nguyen-Widrow method for initializing weights.

Solving the XOR problem using Levenberg-Marquardt

Additional Notes

Some versions of the Levenberg-Marquardt algorithm solve the equation (J^tJ + λ diag(J^tJ) I)δ = J^tE instead of (J^tJ + λI)δ = J^tE, effectively replacing the identity matrix with the diagonal of the approximated Hessian for the weight update rule. According to Wikipedia, this was suggested by Marquardt to incorporate some local curvature estimation. However, for some reason, more recent implementations that use the Levenberg-Marquardt tag do not include Marquardt’s suggestion. Other implementations suggest using an initial λ related to the size of the elements of the approximated Hessian by making λ₀ = t max(diag(J^tJ)), where t is a value chosen by the user.
It is also important to note that the Levenberg-Marquardt is extremely dependent on the initial guess for the network parameters. Depending on the initial weights of the network, the algorithm may converge to a local minima or not converge at all.
Besides that, it is an extremely fast method for neural network learning when compared to the standard backpropagation algorithm.

Finally, if you have any comments about the article or about the code, please let me know it! All suggestions are welcome.

References

Levenberg-Marquardt
- Sam Roweis; Levenberg-Marquardt Optimization, available on: http://www.cs.nyu.edu/~roweis/notes/lm.pdf, 1999.
- Sarle, W.; Neural Network FAQ, part 2 of 7: Learning, available on: http://www.cs.otago.ac.nz/nnweb/FAQ2.html, 2002.
- Wilamowski, B., Chen. Y; Efficient Algorithm for Training Neural Networks with one Hidden Layer, available on: http://cs.olemiss.edu/~ychen/publications/conference/chen_ijcnn99.pdf, 1999.
- Pearlmutter, B. A.; Fast Exact Multiplication by the Hessian, available on: http://eprints.kfupm.edu.sa/40648/1/40648.pdf, 1999.
- Guidry, T. F.; Levenberg Marquardt Using Numerical Derivatives in C#, available in: http://www.trentfguidry.net/post/2009/08/12/Levenberg-Marquardt-using-numerical-derivatives-in-C.aspx, 2009.
Bayesian Regularization
- MacKay, D.; Bayesian methods for neural networks – FAQ, available on: http://www.inference.phy.cam.ac.uk/mackay/Bayes_FAQ.html, 2004.
- Foresee, F. D.; Hagan, M. T.; Gauss-Newton Approximation to Bayesian Learning, available on: http://hagan.ecen.ceat.okstate.edu/icnn97a.pdf, 1997.
- Poland, J.; On the Robustness of Update Strategies for the Bayesian Hyperparameter alpha, available on: http://www-alg.ist.hokudai.ac.jp/~jan/alpha.pdf, 2001.
- The MathWorks, Inc; Neural Network Toolbox – Levenberg-Marquardt (trainlm), available on: http://matlab.izmiran.ru/help/toolbox/nnet/backpr12.html, 2005.
- MacKay, D.; A Practical Bayesian Framework for Backpropagation Networks,
  available on: http://authors.library.caltech.edu/13793/1/MACnc92b.pdf, 1992.

How to become a computer expert!

November 12, 2009 / cesarsouza / Leave a comment

Universal Approximation Theorem

November 5, 2009 / cesarsouza / Leave a comment

The Universal Approximation Theorem states that:

Let φ(·) be a nonconstant, bounded, and monotonically-increasing continuous function. Let I_m₀ denote the m₀-dimensional unit hypercube [0,1]^m₀. The space of continuous functions on I_m₀ is denoted by C(I_m₀). Then, given any function f Э C(I_m₀) and є > 0, there exist an integer m₁ and sets of real constants α_i, b_i and w_ij, where i = 1, …, m₁ and j = 1, …, m₀ such that we may define:

$F( x_1 , dots, x_{m_0} ) = sum_{i=1}^{m_1} alpha_i varphi left( sum_{j=1}^{m_0} w_{i,j} x_j + b_iright)$

as an approximate realization of the function f; that is,

$| F( x_1 , dots, x_{m_0} ) - f ( x_1 , dots, x_{m_0} ) | < varepsilon$

for all x₁, x₂, …, x_m₀ that lie in the input space.

[ From Wikipedia ]

Beautiful, isn’t? 😛

Análise de Componente Principal em C#

November 2, 2009 / cesarsouza / 2 Comments

For the english version, Principal Component Analysis in C#, click this link.

A Análise de Componente Principal (PCA) é uma ferramenta exploratória criada por Karl Pearson em 1901 para identificar tendencias desconhecidas em um conjunto de dados multidimensional. A análise envolve um procedimento matemático que transforma um conjunto de variaveis possivelmente correlatas em um grupo menor de variaveis incorrelatas de denominadas componentes principais. Os primeiros componentes principais correspondem a maior variabilidade possivel contida nos dados, enquanto os componentes posteriores correspondem ao restante da informacao, em ordem decrescente de variabilidade.

Visão Geral

A técnica PCA essencialmente rotaciona o conjunto de pontos ao redor de sua média objetivando seu alinhamento com os primeiros componentes principais. Isto move a maior parte da variancia possivel (usando uma transformação linear) para as primeiras dimensões. Os valores nas dimensoes restantes, portanto, tendem a ser altamente correlacionados e podem ser ignorados com perda mínima de informação.

Para uma explicação mais detalhada sobre PCA, leia o excelente texto Tutorial On Principal Component Analysis, por Lindsay Smith (2002).

AForge.NET Framework

Minha idéia inicial era implementar a PCA internamente no AForge.NET Framework. O AForge.NET é um excelente framework de Inteligência Artificial e Visão Computacional de código aberto escrito para .NET e desenvolvido principalmente por Andrew Kirillov.

No entanto, como isto envolveria muitas mudanças e adições ao framework, acredito que seria muito difícil rever todas mudanças requeridas para adicionar esta funcionalidade diretamente em seu código fonte.

Devido a isso, este novo, e mais limpo, código, foi implementado usando somente os binários do AForge como ponto de partida. Espero que este código ajude outros precisando realizar a Análise de Componente Principal em seus próprios projetos em C#.

Esta nova biblioteca, a que chamei de Accord.NET, extende o AForge.NET adicionando novas features como a Principal Component Analysis, decomposições numéricas e algumas outras poucas transformações matemáticas e ferramentas. Na verdade, esta extensão é muito mais como um campo de teste para novas features que eu gostaria de ver em versões futuras do AForge. Mas para maior simplicidade, omiti as partes não relacionadas com PCA e criei esta versão compacta apenas para suportar a análise.

Decisões de Implementação

Como pessoas que querem utilizar PCA em seus projetos geralmente já tem suas próprias classes de matrizes, decidi evitar usar implementações específicas para tornar o código mais flexível. Também tentei evitar dependencias em outros métodos sempre que possível, para tornar o código bastante independente. Acredito que isto também tenha tornado o código mais simples de se entender.

O código está dividido entre estes dois projetos:

Accord.Math, que fornece as ferramentas matemáticas, como decomposições e transfomações; e
Accord.Statistics, que fornece análises e ferramentas estatísticas.

Ambos dependem do núcleo do AForge.NET. Adicionalmente, sua estrutura e organização interna também tenta similar ao máximo a empregada no AForge sempre que possível.

O código fonte fornecido não inclui o código fonte na íntegra do Accord Framework, que permanece como um framework de testes do que seria interessante ver no AForge, mas inclui somente porções limitadas do código necessárias para suportar PCA.

Visão Geral do Código

Below is the main code behind the PCA.

        <span>/// &lt;summary&gt;Computes the Principal Component Analysis algorithm.&lt;/summary&gt;</span><br />        <span>public</span> <span>void</span> Compute()<br />        {<br />            <span>int</span> rows = sourceMatrix.GetLength(0);<br />            <span>int</span> cols = sourceMatrix.GetLength(1);<br /><br />            <span>// Create a new matrix to work upon</span><br />            <span>double</span>[,] matrix = <span>new</span> <span>double</span>[rows, cols];<br /><br /><br />            <span>// Prepare the data, storing it in the new matrix.</span><br />            <span>if</span> (<span>this</span>.analysisMethod == AnalysisMethod.Correlation)<br />            {<br />                <span>for</span> (<span>int</span> i = 0; i &lt; rows; i++)<br />                    <span>for</span> (<span>int</span> j = 0; j &lt; cols; j++)<br />                        <span>// subtract mean and divide by standard deviation (convert to Z Scores)</span><br />                        matrix[i, j] = (sourceMatrix[i, j] - columnMeans[j]) / columnStdDev[j];<br />            }<br />            <span>else</span><br />            {<br />                <span>for</span> (<span>int</span> i = 0; i &lt; rows; i++)<br />                    <span>for</span> (<span>int</span> j = 0; j &lt; cols; j++)<br />                        <span>// Just center the data around the mean. Will have no effect if the</span><br />                        <span>//  data is already centered (the mean will be zero).</span><br />                        matrix[i, j] = (sourceMatrix[i, j] - columnMeans[j]);<br />            }<br /><br /><br /><br />            <span>// Perform the Singular Value Decomposition (SVD) of the matrix</span><br />            SingularValueDecomposition singularDecomposition = <span>new</span> SingularValueDecomposition(matrix);<br />            singularValues = singularDecomposition.Diagonal;<br /><br /><br />            <span>// Eigen values are the square of the singular values</span><br />            <span>for</span> (<span>int</span> i = 0; i &lt; singularValues.Length; i++)<br />            {<br />                eigenValues[i] = singularValues[i] * singularValues[i];<br />            }<br /><br /><br />            <span>//  The principal components of 'Source' are the eigenvectors of Cov(Source). Thus if we</span><br />            <span>//  calculate the SVD of 'matrix' (which is Source standardized), the columns of matrix V</span><br />            <span>//  (right side of SVD) will be the principal components of Source.                        </span><br /><br /><br />            <span>// The right singular vectors contains the principal components of the data matrix</span><br />            <span>this</span>.eigenVectors = singularDecomposition.RightSingularVectors;<br /><br />            <span>// The left singular vectors contains the scores of the principal components</span><br />            <span>this</span>.resultMatrix = singularDecomposition.LeftSingularVectors;<br /><br /><br />            <span>// Calculate proportions</span><br />            <span>double</span> sum = 0;<br />            <span>for</span> (<span>int</span> i = 0; i &lt; eigenValues.Length; i++)<br />                sum += eigenValues[i];<br />            sum = (1.0 / sum);<br /><br />            <span>for</span> (<span>int</span> i = 0; i &lt; eigenValues.Length; i++)<br />                componentProportions[i] = eigenValues[i] * sum;<br /><br /><br />            <span>// Calculate cumulative proportions</span><br />            <span>this</span>.componentCumulative[0] = <span>this</span>.componentProportions[0];<br />            <span>for</span> (<span>int</span> i = 1; i &lt; <span>this</span>.componentCumulative.Length; i++)<br />            {<br />                <span>this</span>.componentCumulative[i] = <span>this</span>.componentCumulative[i - 1] + <span>this</span>.componentProportions[i];<br />            }<br /><br /><br />            <span>// Creates the object-oriented structure to hold the principal components</span><br />            PrincipalComponent[] components = <span>new</span> PrincipalComponent[singularValues.Length];<br />            <span>for</span> (<span>int</span> i = 0; i &lt; components.Length; i++)<br />            {<br />                components[i] = <span>new</span> PrincipalComponent(<span>this</span>, i);<br />            }<br />            <span>this</span>.componentCollection = <span>new</span> PrincipalComponentCollection(components);

/// <summary>Computes the Principal Component Analysis algorithm.</summary> public void Compute() { int rows = sourceMatrix.GetLength(0); int cols = sourceMatrix.GetLength(1); // Create a new matrix to work upon double[,] matrix = new double[rows, cols]; // Prepare the data, storing it in the new matrix. if (this.analysisMethod == AnalysisMethod.Correlation) { for (int i = 0; i < rows; i++) for (int j = 0; j < cols; j++) // subtract mean and divide by standard deviation (convert to Z Scores) matrix[i, j] = (sourceMatrix[i, j] - columnMeans[j]) / columnStdDev[j]; } else { for (int i = 0; i < rows; i++) for (int j = 0; j < cols; j++) // Just center the data around the mean. Will have no effect if the // data is already centered (the mean will be zero). matrix[i, j] = (sourceMatrix[i, j] - columnMeans[j]); } // Perform the Singular Value Decomposition (SVD) of the matrix SingularValueDecomposition singularDecomposition = new SingularValueDecomposition(matrix); singularValues = singularDecomposition.Diagonal; // Eigen values are the square of the singular values for (int i = 0; i < singularValues.Length; i++) { eigenValues[i] = singularValues[i] * singularValues[i]; } // The principal components of 'Source' are the eigenvectors of Cov(Source). Thus if we // calculate the SVD of 'matrix' (which is Source standardized), the columns of matrix V // (right side of SVD) will be the principal components of Source. // The right singular vectors contains the principal components of the data matrix this.eigenVectors = singularDecomposition.RightSingularVectors; // The left singular vectors contains the scores of the principal components this.resultMatrix = singularDecomposition.LeftSingularVectors; // Calculate proportions double sum = 0; for (int i = 0; i < eigenValues.Length; i++) sum += eigenValues[i]; sum = (1.0 / sum); for (int i = 0; i < eigenValues.Length; i++) componentProportions[i] = eigenValues[i] * sum; // Calculate cumulative proportions this.componentCumulative[0] = this.componentProportions[0]; for (int i = 1; i < this.componentCumulative.Length; i++) { this.componentCumulative[i] = this.componentCumulative[i - 1] + this.componentProportions[i]; } // Creates the object-oriented structure to hold the principal components PrincipalComponent[] components = new PrincipalComponent[singularValues.Length]; for (int i = 0; i < components.Length; i++) { components[i] = new PrincipalComponent(this, i); } this.componentCollection = new PrincipalComponentCollection(components);

Usando o Código

Para realizar uma análise simples, você pode simplesmente instanciar um novo objeto PrincipalComponentAnalysis passando seu conjunto de dados e chamar seu método Compute para computar o modelo. Então você pode simplesmente chamar o método Transform para projetar os dados no espaço dos componentes principais.

Abaixo está um exemplo de código demonstrando seu uso.

            <span>// Creates the Principal Component Analysis of the given source</span><br />            PrincipalComponentAnalysis pca = <span>new</span> PrincipalComponentAnalysis(sourceMatrix,<br />                PrincipalComponentAnalysis.AnalysisMethod.Correlation);<br />    <br /><br />            <span>// Compute the Principal Component Analysis</span><br />            pca.Compute();<br /><br />            <span>// Creates a projection considering 80% of the information</span><br />            <span>double</span>[,] components = pca.Transform(sourceMatrix, 0.8f, <span>true</span>);

// Creates the Principal Component Analysis of the given source PrincipalComponentAnalysis pca = new PrincipalComponentAnalysis(sourceMatrix, PrincipalComponentAnalysis.AnalysisMethod.Correlation); // Compute the Principal Component Analysis pca.Compute(); // Creates a projection considering 80% of the information double[,] components = pca.Transform(sourceMatrix, 0.8f, true);

Demonstração

Para demonstrar o uso da Análise de Componente Principal, criei uma simples aplicação Windows Forms que efetua uma análise de dados simples e as transformações PCA em planilhas Excel.

The application can open Excel workbooks. Here we are loading some random Gaussian data, some random Poisson data, and a linear multiplication of the first variable (thus also being Gaussian).

Simple descriptive analysis of the source data, with a histogram plot of the first variable. We can see it fits a Gaussian distribution with 1.0 mean and 1.0 standard deviation.

Here we perform PCA by using the Correlation method (actually the transformation uses SVD on the standardized data, not the correlation matrix, but the effect is the same). Because the third variable is a linear multiplication of the first, the analysis detected it was totally irrelevant, thus having zero importance level.

Now we can make a projection of the source data using only the first two components.

Nota: Os componentes principais não são únicos porque a decomposição em valor singular não é única. Além do mais, os sinais na matriz de rotação são arbitrários, e portanto podem difererir entre diferentes programas para PCA.

junto com a aplicação de demonstração vem uma planilha Excel contendo vários exemplos de dados. O primeiro exemplo é o mesmo utilizado por Lindsay em seu Tutorial on Principal Component Analysis. Outros exemplos incluem dados Gaussianos, dados não correlacionados e combinações lineares de dados Gaussianos para exemplificar melhor a análise.

Bem, novamente, espero que alguém ache este código útil! 😛

César Souza

Month: November 2009

Como enviar cartas pela internet

Para o território nacional

Para outros países

Neural Network Learning by the Levenberg-Marquardt Algorithm with Bayesian Regularization (part 1)

Downloads

Contents

Overview

Neural Networks

AForge Framework

Network training as a function optimization problem

Levenberg-Marquardt Algorithm

Computing the Jacobian

Approximating the Hessian

Solving the Levenberg-Marquardt equation

General Levenberg-Marquardt Algorithm

Limitations

Neural Network Learning by the Levenberg-Marquardt Algorithm with Bayesian Regularization (part 2)

Downloads

Contents

Bayesian Regularization

Expanded Levenberg-Marquardt Algorithm

Source Code

Using the Code

Sample applications

Additional Notes

References

How to become a computer expert!

Universal Approximation Theorem

Análise de Componente Principal em C#

Visão Geral

AForge.NET Framework

Decisões de Implementação

Visão Geral do Código

Usando o Código

Demonstração