I'm working on a (naïve) algorithm for portfolio optimization using GA. It takes a list of stocks, calculates its expected returns and the covariance between all of them and then it returns the portfolio weights that would produce the highest return of investment given a certain maximum risk the investor is willing to be exposed to.

The expected return of the portfolio is calculated multiplying the weight of the stock (a percentage) by that stock's expected return, while the risk is calculated using the matrix equation w^T * Cov * w, where w is the array of weights and Cov is the covariance matrix.

I'm having trouble improving performance with my algorithm. I've timed the execution of the different parts of the algorithm and found that the crossing is the part that takes the longest (around 5 seconds per iteration).


This is the main class. The portfolio is generated by calling its generate_portfolio() method after instantiating the class with the maximum risk, the stock returns and the covariance matrix. It runs for generations generations unless it fails to improve after 100 generations, in which case it short-cuts the loop and returns. This has been the case every time I've executed it.

class Genetic_Algorithm_Portfolio:
    def __init__(self, max_risk, generations, returns, cov_matrix):
        Creates a genetic algorithm portfolio generator

        - population: list of Candidate objects
        - max_risk: risk constraint the portfolio must meet
        - generations: number of generations the algorithm will run to create a better solution
        self.max_risk = max_risk
        self.generations = generations
        self.returns = returns
        self.cov_matrix = cov_matrix

    def generate_portfolio(self):
        #print("Generating portfolio")
        # Map solutions so that the Candidate objects also have normalized fitness value and cummulative sum
        t0 = time.time()
        random_generator = Random_Portfolio_Generator(5000, self.max_risk, self.returns, self.cov_matrix)
        t1 = time.time()
        print("Time to generate random set: " + str(int(t1 - t0)) + " seconds")
        self.initial_population = random_generator.generate_solutions()[0]
        if len(self.initial_population) == 0:
            return None
        current_generation = self.initial_population
        initial_weight = [0]*len(self.initial_population[0].weights)
        best_solution = Candidate(
            portfolioUtilities.calculate_return(initial_weight, self.returns),
            portfolioUtilities.calculate_std(initial_weight, self.cov_matrix)

        # With this initial population, run 100 generations
        improvements = 100
        for i in range(self.generations):
            t_g0 = time.time()
            if improvements == 0:
                print("No improvements after 100 generations")

            # Next generation will be those selected by the selection method + the 10 best by expected return from the current generation
            t_s0 = time.time()
            next_gen_selected = self.next_generation(current_generation)
            current_generation.sort(key=lambda x: x.expected_return, reverse=True)
            next_gen_selected = next_gen_selected + current_generation[:10]
            t_s1 = time.time()
            print(f"Next generation selection: {str(int(t_s1 - t_s0))} seconds")

            # Pair parents
            t_p0 = time.time()
            parents = self.pair(next_gen_selected)
            t_p1 = time.time()
            print(f"Parents selection: {str(int(t_p1 - t_p0))} seconds")
            # Cross them and create the next gen
            next_gen = []
            t_c0 = time.time()
            for j in range(len(parents)):
                next_gen = next_gen + self.cross(parents[j])

            # Check if you found a better solution
            next_gen.sort(key=lambda x: x.expected_return, reverse=True)
            t_c1 = time.time()
            print(f"Cross: {str(int(t_c1 - t_c0))} seconds")
            if len(next_gen) == 0:

            next_gen_best_solution = next_gen[0]
            if (best_solution.expected_return < next_gen_best_solution.expected_return):
                best_solution = next_gen_best_solution
                improvements = 100
                improvements -= 1
            # Next generation is now current generation
            current_generation = next_gen
            t_g1 = time.time()
            print(f">>> Generation {i} time: {str(int(t_g1 - t_g0))} seconds")
        return best_solution

Random solutions

I generate 5000 random solutions simply checking that all weights are above 0 and the risk is under the max risk passed as parameter:

lass Random_Portfolio_Generator:
    def __init__(self, n, max_risk, returns, cov_matrix):
        self.num_solutions = n
        self.max_risk = max_risk
        self.returns = returns
        self.cov_matrix = cov_matrix

    def generate_solutions(self):
        solutions = []
        best_solution = []
        best_return = float("-inf")
        best_std = 0
        seed = 42
        counter = 0
        attempts = 0
        while counter < self.num_solutions and attempts < 1000:
            seed += 1
            weights = self.generate_random_weights()
            solution = Candidate(
                portfolioUtilities.calculate_return(weights, self.returns),
                portfolioUtilities.calculate_std(weights, self.cov_matrix)
            if (solution.std <= self.max_risk):
                attempts = 0
                counter +=1
                if (solution.expected_return > best_return):
                    best_return = solution.expected_return
                    best_solution = solution
                    best_std = solution.std
                attempts += 1
        return solutions, best_solution, best_return, best_std
    def generate_random_weights(self):
        n = len(self.returns)
        arr = [0] * n
        indexes = list(range(n))
        sum = 0;
        for i in range(n - 1):
            rand = random.randint(0,1000 - sum)
            index = random.choice(indexes)
            sum += rand
            arr[index] = rand
        arr[indexes[0]] = 1000 - sum
        return np.array([x / 1000 for x in arr])


My selection method is encapsulated in the next_generation() method, which uses roulette selection by cumulative sum of the expected return of each portfolio. I also sort the current generation by expected return and add the 10 best chromosomes to the next generation.

    def next_generation(self, current):
        # Select chromosomes to pair
        cum_sums = [o.cum_sum for o in self.initial_population]
        selected = []
        for x in range(len(cum_sums)//2):
            selected.append(self.roulette(cum_sums, np.random.rand()))

            # In case the roulette picks an index that has already been picked
            while len(set(selected)) != len(selected):
                selected[x] = self.roulette(cum_sums, np.random.rand())
        return [self.initial_population[int(selected[x])] for x in range(len(self.initial_population)//2)]

    def roulette(self, cum_sum, chance):
        veriable = list(cum_sum.copy())
        veriable = sorted(veriable)
        return veriable.index(chance)


I pair the next generation by a "next best" approach:

    def pair(self, generation):
        # Pairs selected chromosomes by fitness
        return [generation[i:i + 2] for i in range(0, len(generation), 2)]


To cross the pairs, I take a middle point d = (x1 - x2) / 3 and add it and subtract it from each parent to create 4 new solutions. I found this method in a paper after finding that the traditional crossing methods ended up always creating invalid solutions (more risk than the max tolerated by the user). I filter those that don't fit that criteria neither:

    def cross(self, parents):
        x1 = parents[0].weights
        x2 = parents[1].weights
        d = (x1 - x2) / 3
        x3 = self.mutate(x1 + d)
        x4 = self.mutate(x1 - d)
        x5 = self.mutate(x2 + d)
        x6 = self.mutate(x2 - d)

        return sorted(self.judge_candidates([
                portfolioUtilities.calculate_return(x3, self.returns),
                portfolioUtilities.calculate_std(x3, self.cov_matrix)
                portfolioUtilities.calculate_return(x4, self.returns),
                portfolioUtilities.calculate_std(x4, self.cov_matrix)
                portfolioUtilities.calculate_return(x5, self.returns),
                portfolioUtilities.calculate_std(x5, self.cov_matrix)
                portfolioUtilities.calculate_return(x6, self.returns),
                portfolioUtilities.calculate_std(x6, self.cov_matrix)
        ]), key=lambda x: x.expected_return, reverse=True)[:2]

    def judge_candidates(self, candidates):
        return list(filter(lambda x: x.std <= self.max_risk and all(i >= 0 for i in x.weights), candidates))


Finally, I mutate two random genes in the chromosome with probability 2% in random.gauss(0, 0.001).

    def mutate(self, candidate):
        if (random.randint(0, 100) <= 2):
            gene_one = random.randint(0, len(candidate) - 1)
            gene_two = random.randint(0, len(candidate) - 1)
            deviation = random.gauss(0, 0.001)
            candidate[gene_one] += deviation
            candidate[gene_two] -= deviation
        return candidate


These are the two utility classes I use in this class.

import math

class Candidate:
    def __init__(self, weights, expected_return, std):
        self.weights = weights
        self.expected_return = expected_return
        self.std = math.sqrt(std)
        self.print = f"weights: {self.weights}, return = {self.expected_return}, std = {self.std}"
import numpy as np

def calculate_return(weights, returns):
    return (np.asarray(weights) * returns).sum()

def calculate_std(weights, cov_matrix):
    return np.dot(np.dot(np.asarray(weights), cov_matrix),np.asarray(weights).transpose())

Any ideas or criticisms about what I may be doing wrong?

  • 1
    \$\begingroup\$ @Emma I don't know if I can write LaTeX here, but the problem is: maximize Sum(i=1, n) w_i * r_i subjected to w^T * Cov * w <= max_risk, Sum(i=1, n) w_i = 1, w_i >= 0, where w is the weight vector, r_i is the return for the i-th stock, Cov is the covariance matrix between stocks and max_risk is the maximum risk the client is willing to tolerate. \$\endgroup\$
    – Heathcliff
    Nov 27, 2019 at 1:23
  • 3
    \$\begingroup\$ I barely know any Python, so it won't risk posting an answer, but I have to say this is probably some of the most readable Python I've ever seen. About the algorithm itself though... your population size (num_solutions, 5000) seems way too high. I've found personally that in general, sizes from 20 to 100 work best, or about 10 times your number of dimensions, but I don't know how many dimensions you have here. \$\endgroup\$
    – cliesens
    Nov 27, 2019 at 2:03
  • 2
    \$\begingroup\$ Some other things you could try : use other selection methods (again, personally, I've found roulette selection to not be the most efficient in most cases) such as tournament selection, or some custom ones (keep the best n% for crossover,...). They shouldn't be too hard to implement, and they will often provide the most noticeable changes in performance. Finally, you could also try making your mutation rate vary over time, although I would suggest trying this only after you've noticed actual improvements in your algorithm. \$\endgroup\$
    – cliesens
    Nov 27, 2019 at 2:10
  • \$\begingroup\$ @cliesens You are not strictly required to know the programming language as long as you can give feedback with regard to the algorithm. \$\endgroup\$
    – AlexV
    Nov 27, 2019 at 8:02
  • \$\begingroup\$ The current question title, which states your concerns about the code, is too general to be useful here. Please edit to the site standard, which is for the title to simply state the task accomplished by the code. Please see How to get the best value out of Code Review: Asking Questions for guidance on writing good question titles. \$\endgroup\$ Nov 27, 2019 at 11:28

1 Answer 1


Few tips:

  • try hashing fitness function, create dictionary of hashed individuals and their fitness values; fitness evaluation is often costly; no reason to evaluate same specimen twice
  • don't use roulette selection, you can loose selective pressure, when there's little difference between fitnesses or converge too fast, when there's huge deviation; use tournament selection instead, it doesn't care about fitness distribution and maintains selective pressure well
  • GA is rather well suited for parallel computation; exploit that fact; you got nice multiprocessing package in python's standard lib
  • try to use hill climber to improve some of the individuals; there was a lot of research about hybrids: simulated annealing, tabu search, FIHC. It depends on your problem and encoding though.

Now most important stuff - vanilla GA is slow, because it searches through space like a blind man. Try to learn about linkage learning (discrete problems) and fitness landscape analysis (continous problems). You'll be able to crash your problem way much faster. If you'd like to know more about it, get some papers, please create separate thread on https://datascience.stackexchange.com/ or https://stats.stackexchange.com/. We can talk there.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.