Friday, April 6, 2012

Generate random matrices with a given condition number (or spectrum)

It is really easy to perform this task. Just generate a random matrix, for example, in Matlab:
n=4
A=rand(n)
then we need no extract the (random) basis with the QR decomposition.
[Q,R]=qr(A)
and then we fix the spectrum. The best possibly conditioned matrix would have a condition number of one. Therefore we can achieve that with
D=diag(ones(n,1))
M=Q*D*Q'
This would produce a normal matrix $MM^T = M^T M$. If we want a non-normal matrix, we just need another random basis, that we can take from another random matrix $B$.
B=rand(n)
[V,R]=qr(A)
M=Q*D*V'
If we want to fix only the first (largest) eigenvalue to one, but leaving a random spectrum, we generate a random sequence
sp=randn(n,1)
sp=sp/max(sp)
D=diag(sp)
M=Q*D*V'

Sunday, April 1, 2012

Practical kernel trick

When using an SVM, a naive usage of the software may lead us to use the standard kernels that are offered to us and perform the analysis as well as we can.

However, one can choose to build his own kernels given some prior information about the problem. For example, we could have a dataset with values of instantiated variables for each datum, but we might also have correlations telling us when those data simultaneously appear. For example, in credit rating, one may have variables such as income for each of the customers. One might also have some kind of information linking customers, such as marital status, housing status, professional category and so on. Therefore, apart from the points in $\mathbb{R}^n$ (wages), we also have some kind of coexpression graph, where two individuals score higher with each other if their categorical variables coincide, so that we build a matrix $S_{i,j}=\mbox{#}(k,j)$, where $\mbox{#}(i,j)$ is the coincidence of the vales for all variables between individual $i$ and individual $j$. We then could make a custom kernels in this way:
$$
K=K_{\sigma} + \gamma S^*
$$
Where $K_{\sigma}$ is the classical Gaussian kernel build from the euclidean distances of the wages, with parameter $\sigma$, and $S^*$ is a projection of the previous $S$ onto the cone of positive semi-definite matrices. $\gamma$ is a tuning parameter. Now, how can we use it with our SVM software?

LibSVM can be fed with the previous kernel matrix $K$ if you work with the C/C++ or the Matlab interface. Just make sure to insert an index from one to the number of data as the first column. If you are working with R with the library e1071, you won't be able to do this.

It is well known to practitioners that the kernel-trick can help us in this case. Since the Mercer's kernel expansion is
$$
k(x,y)=\sum_i \lambda_i \phi_i(x) \phi_i(y) = \Phi(x)^T \Phi(y)
$$
Where $\lambda_i$ and $\phi_i$ are the eigenvalues and eigenfunctions of the integral operator whose positive definite kernel is $k$. $\Phi(x) = [\sqrt{\lambda_1} \phi_1(x), \sqrt{\lambda_2} \phi_2(x),\ldots]$ is, therefore, the "lifting" map from $\mathbb{R}^n$ to the Hilbert space of square-summable infinite sequences $l^2$. Since we are working with a finite approximation to this operator, the eigen-approximation is given by the eigen values and eigenvectors of the kernel matrices above. We will always be able to compute these. Therefore, we substitute $ \Phi(x)^T \Phi(y)$ by its finite approximation using the eigenvalues and eigenvectors of the matrix, so that using a standard linear kernel with the lifting map (the images under the approximated $\Phi$) we achieve the values in $K$.

The following example is in R. Assume we want to compute some credit score input, and the individuals' salaries are the ones given in the variable salaries. The experience with these customers is +1 if they were good customers, and -1 if they were bad customers, and stored in the variable y.
library(e1071)
salaries=c(20000,30000,55000,50000,30000,25000)
y=c(+1,+1,+1,-1,-1,-1)
 It is clear that the experience and the salaries are not very correlated. Imagine betting is a huge problem for some segment of the population, and we have the individuals' record of using betting facilities, in this matrix
betting=c(0,0,0,1,1,1)

S=betting%*%t(betting)
Seig=eigen(S)
Seig$vectors
values=Seig$values
values[values<0]=0
D=diag(values)
Ss=Ss=U%*%D%*%t(U)
s=0.000000005
Ks=as.matrix(exp(-s*dist(salaries)^2))
diag(Ks)<-1
The gaussian distances are:  
2 0.606530660                                               
3 0.002187491 0.043936934                                   
4 0.011108997 0.135335283 0.882496903                       
5 0.606530660 1.000000000 0.043936934 0.135335283           
6 0.882496903 0.882496903 0.011108997 0.043936934 0.88249690


We now see that the true correlations between these individuals are not very strong. Therefore, supplementing Ks with S is a good idea. Ks+gamma*S is:
            1          2           3          4          5          6
1 0.000000000 0.60653066 0.002187491 0.01110900 0.60653066 0.88249690
2 0.606530660 0.00000000 0.043936934 0.13533528 1.00000000 0.88249690
3 0.002187491 0.04393693 0.000000000 0.88249690 0.04393693 0.01110900
4 0.011108997 0.13533528 0.882496903 1.00000000 1.13533528 1.04393693
5 0.606530660 1.00000000 0.043936934 1.13533528 1.00000000 1.88249690
6 0.882496903 0.88249690 0.011108997 1.04393693 1.88249690 1.00000000


We see that there is some more correlation between the betting public. Following the R program
K=Ks+gamma*S
Keig=eigen(K)
phi=Keig$vectors%*%diag(sqrt(Keig$values))
Now we train the SVM with this embedding, or lifting map, $\Phi$. Observe that if $K=VEV^T=VE_q E_q V^T = \Phi^T \Phi$ then the embedding is $\Phi^T=VE_q$, where $E_q$ has the square root of the eigenvalues.
svm(phi, y, type="C", kernel="linear", cost=1)
Subsequent customers (after retraining) will be given a more fitting credit rating in this fashion than using only their salaries.

This trick is also the basis of the celebrated technique of random features (check out the paper).

Saturday, March 31, 2012

Scientific computing at home with CUDA

I wanted to write about CUDA because I feel it is the future of non-desktop computation, meaning any process that goes beyond data movement and has to apply a scientific computation to obtain information from raw data.

CUDA (Compute Unified Device Architecture) is a form of GPGPU that uses the graphics processors on NVIDIA graphics cards as computation units. One can send programs (called kernels) that perform those computations to the graphics device.

The importance of CUDA lies in the specific hardware architecture aimed at performing vector computations. Scientific and graphics software making extensive use of arithmetic operations will therefore benefit from CUDA parallelization (this includes everywhere you see matrix algebra, such as in quadratic optimization including SVMs, PCA, ICA, CCA, and other discretized operations such as fast Fourier transform, wavelet filter banks and so on). On the other hand, software using large memory transfers are impervious to (any kind of arithmetic) parallelization (this includes databases, web servers, algorithmic/protocol-driven networking software, etc.).

It is interesting for data processing practitioners in the sense that it can cope with the larga datasets. Modern CPUs contain a high among of L2 cache, and the tree-structured cache may expand up to three levels. This design is so because software programs have high data coherency, meaning predominance of serial (and thus continuously accessible) code. However, when working quickly across a large dataset, this cache design performance decays very rapidly. The GPU memory interface is very different from that of the CPU. GPUs use massive parallel interfaces in order to connect with their memory. For example, the GTX 280 uses a 512-bit interface to its high performance GDDR3 memory. This interface is approximately 10 times faster than a typical CPU-RAM interface. On the other hand, GPUs lack the vast amount of memory that the main CPU system enjoy. Customized, large memory parallel GPUs are commercially available at a higher price than that of a home high-performing gaming system. But a nice approach to higly-parallelized software can be taken notwithstanding.

There is a serious drawback from a software engineering point of view, and that is that, Unlike most programming languages, CUDA is coupled very closely together with the hardware implementation. While CPU families do not basically change to maintain retro-compatibility, CUDA-enabled hardware significantly change in basic architectural design.

An important concept to be aware of is the thread hierarchy. CPUs are designed to run just a few threads. GPUs, on the other hand, are designed to process thousands of threads simultaneously. So, in order to take full advantage of your graphics card, the problem must be liable to be broken down to small pieces. Important CUDA tree-like thread organizational structures one has to bear in mind:
Half-Warp: A half-warp is a group of 16 consecutive threads. Half-warp threads are generally executed together and aligned, for instance, threads 0-15 will be in the same half-warp, 16-31, and so on.
Warp: A warp of threads is a group of 32 consecutive threads. On future computing devices from NVIDIA, it might be possible that all threads in the same Warp are generally executed together in parallel. Therefore, it is a good idea to make your programs as if all threads within the same warp will execute together in parallel.
Block: A block is a collection of threads. For technical reasons, blocks should have at least 192 threads to obtain maximum efficiency and full latency avoidance. Typically, blocks might contain 256 threads, 512 threads, or even 768 threads. Here’s the important thing you need to know. Threads within the same block can synchronize with each other, and quickly communicate with each other.
Grid: A grid is a collection of blocks. Blocks can not synchronize with each other, and therefore threads within one block can not synchronize with threads in another block.

Regarding memory, CUDA uses the following kinds:

Global Memory: Global memory can be thought of as the physical memory on your graphics card. All threads can read and write to Global memory.
Shared Memory: A GPU consists of many processors, or multiprocessors. Each multiprocessor has a small amount of shared memory, with a usual size of 16KB. Shared memory is generally used as a very quick working space for threads within a block. It is allocated on a block by block basis. For example, you may have three blocks running consecutively on the same multiprocessor. This means that the maximum amount of shared memory the blocks can reserve is 16KB/3. Threads within the same block can quickly and easily communicate with each other by writing and reading to the shared memory. It’s worth mentioning that the shared memory is at least 100 times faster than global memory, so it’s very advantageous if you can use it correctly.
Texture Memory: A GPU also has texture units and memory which can be taken advantage of in some circumstances. Unlike global memory, texture memory is cached, and is generally read only. This can be taken advantage by using this memory to store chunks of raw data to process.

From the software development point of view, the programs run within a thread are:

Kernels: In a Warp, a kernel is executed multiple times in a SIMD (single instruction multiple data) fashion, which means that the flow of execution in the processors is the same, executing the same instruction, but each processor operation that instruction on different data. We are interested in splitting the software in pieces so that a piece is executed multiple types on different data. For example, in the matrix multiplication A*B, a single kernel may have the task to multiply a row from A and a column from B. When instantiating this kernel in N x M threads, each instantiation (thread) performs the dot-product of a specific row of A and a specific row of B.

Since, at the time of compilation and linkage, the kernel, which is CUDA code, is stored as data in the program (as seen by the operating system), and is subsequently sent to the GPU via the corresponding DMA channel. Once in the GPU, the kernel is instantiated on each processor.

I made a quick Example that can directly be compiled with Visual C++ with help of the the included project files (assuming one has the CUDA toolkit, driver, and Visual C++ environment setup). It can also be compiled anywhere one has a working C++ and NVCC environment.

I will write a follow up article with the basic setip with Visual C++ 2008 Express, which requires some tweaks to get up and running.

Tuesday, March 20, 2012

Multiclass SVM with e1071

When dealing with multi-class classification using the package e1071 for R, which encapsulates LibSVM, one faces the problem of correctly predicting values, since the predict function doesn't seem to deal effectively with this case. In fact, testing the very example that comes in the svm help (?svm on the R command line), one sees the failing performance of the function (albeit working with a correctly fitted model).

data(iris)
attach(iris)

model <- svm(Species ~ ., data = iris)

x <- subset(iris, select = -Species)
y <- Species
model <- svm(x, y, probability = TRUE)

This model is correctly fitted. However

pred <- predict(model, x)
Does not show the correct values

> pred     1      2      3      4      5      6      7      8      9     10     11     12     13     14     15 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    16     17     18     19     20     21     22     23     24     25     26     27     28     29     30 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    31     32     33     34     35     36     37     38     39     40     41     42     43     44     45 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    46     47     48     49     50     51     52     53     54     55     56     57     58     59     60 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    61     62     63     64     65     66     67     68     69     70     71     72     73     74     75 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    76     77     78     79     80     81     82     83     84     85     86     87     88     89     90 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
    91     92     93     94     95     96     97     98     99    100    101    102    103    104    105 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
   106    107    108    109    110    111    112    113    114    115    116    117    118    119    120 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
   121    122    123    124    125    126    127    128    129    130    131    132    133    134    135 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
   136    137    138    139    140    141    142    143    144    145    146    147    148    149    150 
setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa setosa 
Levels: setosa versicolor virginica
We see that every prediction is setosa, even though the dataset is equally divided between the three classes (setosa, versicolor and virginica). It seems that something went wrong with the direct prediction, but one way to overcome this problem is to use the predicted probabilities, that seem to be well computed:

pred <- predict(model, x, decision.values = TRUE, probability = TRUE)
Observe how they look

> attr(pred, "probabilities")         setosa  versicolor   virginica
1   0.980325653 0.011291686 0.008382661
2   0.972927977 0.018061300 0.009010723
3   0.979044327 0.011921880 0.009033793
...
48  0.977140826 0.013050710 0.009808464
49  0.977831001 0.013359834 0.008809164
50  0.980099521 0.011501036 0.008399444
51  0.017740468 0.954734399 0.027525133
52  0.010394167 0.973918376 0.015687457
...
97  0.009263806 0.986123276 0.004612918
98  0.008503724 0.988168405 0.003327871
99  0.025068812 0.965415124 0.009516064
100 0.007514580 0.987584706 0.004900714
101 0.012482541 0.002502134 0.985015325
...
149 0.013669944 0.017618659 0.968711397
150 0.010205071 0.140882630 0.848912299
Now we see that the most probable class is indeed the ground truth and we can correctly classify with the following function
predsvm<-function(model,newdata)
{
  prob<-attr(predict(model, newdata, probability = TRUE),"probabilities")
  n<-dim(prob)[1]
  m<-dim(prob)[2]
 
  me<-which(prob==apply(prob,1,max))
  return(as.factor(model$labels[floor((me-1)/n)+1]))
}
One might also program the following function, that deals with the way the support vector coefficients are stored in the model object, in model$coefs and model$rho:
## Linear Kernel function
K <- function(i,j) crossprod(i,j)

predsvm <- function(object, newdata) {
  ## compute start-index
  start <- c(1, cumsum(object$nSV)+1)
  start <- start[-length(start)]

  ## compute kernel values
  kernel <- sapply (1:object$tot.nSV,
                    function (x) K(object$SV[x,], newdata))

  ## compute raw prediction for classifier (i,j)
  predone <- function (i,j) {
    ## ranges for class i and j:
    ri <- start[i] : (start[i] + object$nSV[i] - 1)
    rj <- start[j] : (start[j] + object$nSV[j] - 1)
    
    ## coefs for (i,j):
    coef1 <- object$coefs[ri, j-1]
    coef2 <- object$coefs[rj, i]

    ## return raw values:
    crossprod(coef1, kernel[ri]) + crossprod(coef2, kernel[rj])
  }

  ## compute votes for all classifiers
  votes <- rep(0,object$nclasses)
  c <- 0 # rho counter
  for (i in 1 : (object$nclasses - 1))
    for (j in (i + 1) : object$nclasses)
      if (predone(i,j) > object$rho[c <- c + 1])
        votes[i] <- votes[i] + 1
      else
        votes[j] <- votes[j] + 1

  ## return winner (index with max. votes)
  object$levels[which(votes %in% max(votes))[1]]
}

Sunday, March 18, 2012

Smoothing splines and interest rate curves

Yield curves are important in Economics and used by finance professionals to analyze bonds and look for trading opportunities and by economists, to try to understand economic conditions.

The yield is the amount in cash that returns to the owners of a security, for example, if a sovereign bond with maturity at time $T$ is bought at time $t$ for a price $P(t,T)$, the yield would be high if $P(t,T)$ is much less than the maturity price $P(T,T)$ (which would obviously be the risk-free amount received by the lender). A higher yield allows the owner to make more profit from the investment, but it may also mean that the security is being sold at a discount as a result of some risk related to the security. For a recent example of this, Greek bonds, which were sold at such discounts because the market anticipated a default, this made yields soar, with risk-averting investors unloading bonds at the lower prices that risk-taking investors were comfortable buying in. Ultimately everybody took a haircut.

For fixed income mathematics, a set of axioms must be assumed, namely:
  1. The market trades continuously over its trading horizon: it extends from the current time to some distant future time such that the maturities of all the instruments to be valued fall between now and the trading horizon..
  2. The market is efficient: information is available to all traders simultaneously, and every trader makes use of all the available information, also related to the axiom that
  3. There are no arbitrage opportunities: the price of a portfolio is the sum of its constituent parts.
  4. The market is complete: any desired cash flow can be obtained from a suitable self-financing strategy based on a portfolio of discount bonds.

It is assumed further that there are no legal barriers to trading and that the market is rational (traders try to maximize their profits).

Now, define the zero yield in terms of the bond price above. The zero yield, as seen at time $t$, of a bond that matures at time $T$, is denoted by $y(t,T)$ and is defined by
$$
P(t,T)=e^{-(T-t) y(t,T) }
$$
for every $t<T$, which means that the current price of the bond is equal to discounted price at maturity (using compound interest).

Suppose that at time $t$ we enter into a forward contract to deliver at time $T_1$ a bond that will mature at time $T_2$ (obviously $T_2>T_1$). Let the forward price of the bond be denoted by $P(t,T_1,T_2)$. At the same time, a bond that matures at time $T_1$ is purchased at a price of $P(t,T_1)$. Again at time $t$, a bond that matures at time $T_2$ is bought with a price of $P(t,T_2)$. Notice that the axiom specifying that there are no arbitrage opportunities implies that the price of the bond maturing at time $T_2$ must be equal to the product of the price of the bond maturing at time $T_1$ and the forward price $P(t,T_1) P(t,T_1,T_2)$. Let the implied forward rate $f(t,T_1,T_2)$, valued at time $t$, for the period $(T_1,T_2)$ as
$$
P(t,T_1,T_2)=e^{-(T_2-T_1) f(t,T_1,T_2) }
$$
Notice that $P(t,T_2) =P(t,T_1) P(t,T_1,T_2)=e^{-(T_1-t) y(t,T_1) } e^{-(T_2-T_1) f(t,T_1,T_2) }$ and, summing the exponents, we get
$$
f(t,T_1,T_2)= \frac{(T_2-t) y(t,T_2) - (T_1 - t) y(t,T_1)}{T_2 - T_1}
$$
which is the period forward rate. However, the instantaneous forward rate is of much greater importance in the theory of the term structure. The instantaneous forward rate for time $T$, as seen at time $t$, is denoted by $f(t,T)$ and is the continuously compounded rate defined by
$$
f(t,T)= \lim_{h \rightarrow 0} {f(t,T,T+h) } = y(t,T) + (T-t) y_T (t,T)
$$
where $y_T = \frac{\partial}{\partial T} y(t,T)$. We can interpret the previous equation as the current  yield value plus the instantaneous change of the yield with the maturity time.

To interpolate curves, and thus have values for all the points within the interval and not only the data points that are made available to us by the problem, we can use a number of different methods. One needs some complexity to be able to capture what the nature is saying, but at the same time this complexity might be cause by some "noise", or perturbations that we want to avoid (market panics, manipulations or poorly registered data). We can choose a parametric family, such as polynomials, fix the order and fit the coefficients. However, this appears as somehow arbitrary. One can choose to interpolate with a non-parametric method such as splines. When using some functional space, one must restrict himself to functions that meet some criteria, not only that approximate well (or maybe the best approximation in that space). This is to avoid overfitting the said noise.

Now we can choose to fit an approximating yield or the forward curve that approximates the true yield or forward curves, with the information given by actual bond prices as of time $t$, by minimizing the following functional (assuming we fit the yield curve $\varphi$):
$$
R(\varphi,t) =  \sum_{i=1}^n (P(t,T_i) - P(t,T_i,\varphi))^2 + \lambda \int \left( \frac{\partial^2}{\partial s^2} \varphi(s) \right)^2 ds.
$$
Where $\lambda$ is a regularization parameter, and $P(t,T_i,\varphi)$ is a functional that prices a bond of maturity $T_i$ at time $t$ with the yield curve $\varphi$. We change formulation, for practical derivation purposes (we now compute the observed yield $y(t,T_i)$ from the observed bond price)
$$
R(\varphi,t) = \sum_{i=1}^n (y(t,T_i) - \varphi(T_i))^2 + \lambda \int \left( \frac{\partial^2}{\partial s^2} \varphi(s) \right)^2 ds.
$$
And then we apply the standard smoothing splines theory to compute the solution, which is:
$$
\mathbf{w} = (I+\lambda K)^{-1} \mathbf{y}
$$
Where $\mathbf{y}$ is a vector with elements $y(t,T_i)$, $K$ is a matrix whose elements are of the form $\int \frac{\partial^2}{\partial s^2} \phi_i (s)  \frac{\partial^2}{\partial s^2} \phi_j (s) ds$, with $\phi$ is a choosen spline basis (notice the dot-product structure) and $\mathbf{w}$ is a vector with the coefficients on that basis, so that the estimated function is
$$
\varphi = \sum_i w_i \phi_i
$$



Thursday, March 15, 2012

Beautiful spectacle in the sky

I enjoy Astronomy so much that I examine the night sky every time I can. Tonight I couldn't wait to describe the delightful sight in the sky I was witnessing: Jupiter and Venus almost at the closest distance (as they meet the eye, not actual Euclidean distance, they were closest two days ago). And these are some beautiful pictures on the Internet (you may see them in space.com) so that the blog records this spectacular celestial art.

 Mostar in Bosnia-Herzegovina

Eagle Harbor, Mich. USA

As per the fabulous site http://www.solarsystemscope.com, the current planet locations are like this, observe the narrow angle in which we can see Jupiter and Venus.

An observer from California at the time of this post would see Venus on the top right of Jupiter in the sunset (in the previous picture of the solar system, put your head horizontally to the left, and see that this would indeed be true).

If the doors of perception were cleansed every thing would appear to man as it is, infinite
William Blake

One feels blessed to having been granted the opportunity to grasp such a marvelous look into infinity.

The coin denomination problem

I read sometime that this is a common question in interviews with Google. The problem states:

Given any set of coins, build an algorithm that returns the minimum amount of coins necessary to make up this number.

The first naive strategy one tends to follow is to use the maximim amount of the maximum valued coins in the set, and then follow with lower valued coins. This usually works with real world coin systems, with follow a pattern that makes this possible. However, the general problem is slightly more complicated. Consider the following set of coins: {1,2,5,8,10}. This set won't work with this strategy.

A better approach consists on considering the lower-valued coins not only when more higher valued coins are substracted, but on every iteration. When we are finished, we perform the same operation recursively, without using the highest-valued coin at the present recurseve call. In the example, we would take the 10-coin out, in the next recursive call, we take the 8-coin out, an so forth. Thus, we take into account all possibilities. Of course, we should track the minimum number of coins to cut the tree.

The following program implements the solution to this problem:

/* This program solves the coin denomination problem.
//
// You may redistribute this file according to the GPL terms.
// Copyright 2006 Javier Arriero Pais
*/

#include <stdio.h>
#include <stdlib.h>

#define MAXC 5
int coins[MAXC] = {1,2,5,8,10};

int amount = 13;

int numcoins (int, int);

int main (int argc, char *argv[])
{
 if (argc &gt; 1)
amount = strtol (argv[1], NULL, 10);

 printf ("Minimal number of coins: %i\n", numcoins(amount, MAXC) );

 return 0;
}


int numcoins (int amount, int maxc)
{
 int i;
 int auxamount=amount, mon=10000, auxmon=0, auxmon2;

 if (maxc&lt;=0) return 0;
 if (maxc==1) return amount;
 if (amount &lt;= 0) return 0;
 if (amount == 1) return 1;

 i=maxc-1;
 while (auxamount &gt;= coins[i]) {
auxamount-=coins[i];
auxmon2 = numcoins (auxamount, maxc-1);
auxmon++;
if (auxmon + auxmon2 &lt; mon) mon = auxmon2 + auxmon;
 }

 auxmon2 = numcoins(amount, maxc-1);

 if (auxmon2 != 0 && auxmon2 < mon)
mon = auxmon2;

 return mon;
}