Arquivo de classificacao - Página 2 de 2

Regressão de Crista do Kernel

19/12/202518/10/2025 Por antonino

Continuando nossa jornada pelo guia do scikit-learn, chegamos a uma técnica poderosa que combina regularização com métodos de kernel: a Regressão de Crista do Kernel. Primordialmente, esta abordagem estende os conceitos de regularização que discutimos anteriormente, incorporando a flexibilidade dos métodos de kernel para lidar com relações não-lineares nos dados.

Conceitos Fundamentais

Analogamente à Regressão de Crista tradicional que exploramos em seções anteriores, a versão com kernel mantém o princípio de regularização, mas opera em um espaço de características transformado através de funções kernel. Similarmente aos métodos de kernel que encontramos em SVM, esta técnica permite capturar padrões complexos sem explicitamente computar transformações de alta dimensionalidade.

Formulação Matemática

Problema de Otimização

A Regressão de Crista do Kernel resolve o seguinte problema de otimização:

\(\min_{\alpha} \|K\alpha – y\|^2_2 + \lambda \alpha^T K\alpha\)

Onde:

K é a matriz do kernel
α são os coeficientes no espaço do kernel
y é o vetor de targets
λ é o parâmetro de regularização

Solução Fechada

A solução pode ser expressa em forma fechada:

\(\alpha = (K + \lambda I)^{-1}y\)

As previsões para novos pontos são dadas por:

\(\hat{y} = K_{test}\alpha\)

Vantagens da Abordagem com Kernel

Conforme demonstramos em técnicas anteriores, a incorporação de kernels oferece benefícios significativos:

Capacidade de modelar relações não-lineares complexas
Operação implícita em espaços de alta dimensionalidade
Flexibilidade através da escolha do kernel
Preservação da estrutura de regularização

Kernels Suportados

O scikit-learn oferece diversos kernels pré-definidos:

linear: \(K(x, x’) = x^T x’\)
polynomial: \(K(x, x’) = (\gamma x^T x’ + r)^d\)
rbf: \(K(x, x’) = \exp(-\gamma \|x – x’\|^2)\)
sigmoid: \(K(x, x’) = \tanh(\gamma x^T x’ + r)\)

Parâmetros e Hiperparâmetros

Parâmetro de Regularização (alpha)

Controla o trade-off entre fitting dos dados e simplicidade do modelo. Analogamente ao parâmetro lambda na regressão ridge tradicional, valores maiores aumentam a regularização.

Parâmetros do Kernel

Cada kernel possui seus próprios parâmetros específicos:

gamma: Controla a influência de cada amostra (RBF, polynomial)
degree: Grau do polinômio (polynomial kernel)
coef0: Termo independente (polynomial, sigmoid)

Conexões com Tópicos Anteriores

Similarmente aos conceitos que exploramos em regressão regularizada e métodos de kernel:

Combina a regularização L2 da Regressão Ridge
Utiliza o truque do kernel como em SVM
Mantém soluções em forma fechada como na regressão linear
Oferece flexibilidade comparável a métodos baseados em árvores

Aplicações Práticas

Primordialmente, a Regressão de Crista do Kernel é útil em cenários como:

Problemas de regressão com relações não-lineares
Datasets com dimensionalidade moderada a alta
Cenários onde interpretabilidade não é crítica
Aplicações que beneficiam de modelos suaves

Exemplo Prático em Python

Para ilustrar a aplicação da Regressão de Crista do Kernel, implementemos um estudo comparativo detalhado:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.kernel_ridge import KernelRidge
from sklearn.linear_model import Ridge
from sklearn.model_selection import train_test_split, GridSearchCV
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import make_regression, make_friedman1
import warnings
warnings.filterwarnings('ignore')
 
'''
Estudo comparativo da Regressão de Crista do Kernel
em diferentes cenários de dados
'''
 
print("=== REGRESSÃO DE CRISTA DO KERNEL ===")
 
# Criando diferentes cenários de dados
print("\n1. CONFIGURAÇÃO DOS CENÁRIOS")
 
# Cenário 1: Relação linear com ruído
X_linear, y_linear = make_regression(
    n_samples=200, n_features=2, noise=10, random_state=42
)
 
# Cenário 2: Relação não-linear (Friedman1)
X_nonlinear, y_nonlinear = make_friedman1(
    n_samples=200, n_features=5, noise=0.1, random_state=42
)
 
# Cenário 3: Dados com interações complexas
np.random.seed(42)
X_complex = np.random.randn(200, 3)
y_complex = (X_complex[:, 0]**2 + 
             np.sin(X_complex[:, 1] * 2) + 
             X_complex[:, 0] * X_complex[:, 2] + 
             np.random.normal(0, 0.1, 200))
 
cenarios = [
    ('Linear', X_linear, y_linear),
    ('Não-Linear (Friedman1)', X_nonlinear, y_nonlinear),
    ('Complexo', X_complex, y_complex)
]
 
'''
Configuração dos modelos comparativos
'''
print("\n2. CONFIGURAÇÃO DOS MODELOS")
 
models = {
    'Ridge Linear': Ridge(),
    'KernelRidge RBF': KernelRidge(kernel='rbf'),
    'KernelRidge Polinomial': KernelRidge(kernel='polynomial'),
    'KernelRidge Linear': KernelRidge(kernel='linear')
}
 
# Espaço de parâmetros para tuning
param_grids = {
    'Ridge Linear': {'alpha': [0.1, 1.0, 10.0, 100.0]},
    'KernelRidge RBF': {
        'alpha': [0.1, 1.0, 10.0],
        'gamma': [0.1, 1.0, 10.0]
    },
    'KernelRidge Polinomial': {
        'alpha': [0.1, 1.0, 10.0],
        'degree': [2, 3, 4],
        'gamma': [0.1, 1.0]
    },
    'KernelRidge Linear': {
        'alpha': [0.1, 1.0, 10.0, 100.0]
    }
}
 
'''
Avaliação sistemática em todos os cenários
'''
print("\n3. AVALIAÇÃO COMPARATIVA")
 
resultados_completos = {}
 
for cenario_nome, X, y in cenarios:
    print(f"\n--- {cenario_nome} ---")
    print(f"Dimensões: {X.shape}")
    
    # Padronizando os dados
    scaler = StandardScaler()
    X_scaled = scaler.fit_transform(X)
    y_scaled = StandardScaler().fit_transform(y.reshape(-1, 1)).ravel()
    
    # Dividindo em treino e teste
    X_train, X_test, y_train, y_test = train_test_split(
        X_scaled, y_scaled, test_size=0.3, random_state=42
    )
    
    cenario_resultados = {}
    
    for model_name in models.keys():
        print(f"  Treinando {model_name}...")
        
        try:
            # Grid Search com validação cruzada
            grid_search = GridSearchCV(
                models[model_name], 
                param_grids[model_name], 
                cv=5, 
                scoring='neg_mean_squared_error',
                n_jobs=-1
            )
            
            grid_search.fit(X_train, y_train)
            
            # Melhor modelo
            best_model = grid_search.best_estimator_
            best_params = grid_search.best_params_
            
            # Previsões e métricas
            y_pred_train = best_model.predict(X_train)
            y_pred_test = best_model.predict(X_test)
            
            mse_train = mean_squared_error(y_train, y_pred_train)
            mse_test = mean_squared_error(y_test, y_pred_test)
            r2_train = r2_score(y_train, y_pred_train)
            r2_test = r2_score(y_test, y_pred_test)
            
            # Armazenando resultados
            cenario_resultados[model_name] = {
                'model': best_model,
                'best_params': best_params,
                'mse_train': mse_train,
                'mse_test': mse_test,
                'r2_train': r2_train,
                'r2_test': r2_test,
                'best_score': -grid_search.best_score_  # Convertendo para MSE positivo
            }
            
            print(f"    Melhores parâmetros: {best_params}")
            print(f"    MSE teste: {mse_test:.4f}, R² teste: {r2_test:.4f}")
            
        except Exception as e:
            print(f"    ERRO: {e}")
            cenario_resultados[model_name] = {'error': str(e)}
    
    resultados_completos[cenario_nome] = cenario_resultados
 
'''
Visualização comparativa dos resultados - VERSÃO CORRIGIDA
'''
print("\n4. VISUALIZAÇÃO DOS RESULTADOS")

# Criando visualizações para cada cenário
for cenario_idx, (cenario_nome, X, y) in enumerate(cenarios):
    # CORREÇÃO: Usar o scaler específico para cada cenário
    scaler_cenario = StandardScaler()
    X_scaled_cenario = scaler_cenario.fit_transform(X)
    y_scaled_cenario = StandardScaler().fit_transform(y.reshape(-1, 1)).ravel()
    
    # Para visualização, usar apenas os dois primeiros features se houver mais
    if X.shape[1] > 2:
        X_viz = X[:, :2]
        print(f"  {cenario_nome}: Usando 2 features para visualização de {X.shape[1]} features totais")
    else:
        X_viz = X
    
    # CORREÇÃO: Ajustar layout dos subplots baseado no número de modelos
    n_models = min(4, len(models))  # Máximo 4 modelos para visualização
    n_rows = 2
    n_cols = 3  # 1 para dados originais + 4 para modelos + 1 para comparação
    
    fig, axes = plt.subplots(n_rows, n_cols, figsize=(18, 10))
    fig.suptitle(f'Regressão de Crista do Kernel - {cenario_nome}', fontsize=16)
    
    # Aplanar axes para facilitar acesso
    if n_rows > 1 and n_cols > 1:
        axes_flat = axes.flatten()
    else:
        axes_flat = [axes] if n_cols == 1 else axes
    
    # Plot 1: Dados originais
    ax_idx = 0
    if X_viz.shape[1] == 1:
        axes_flat[ax_idx].scatter(X_viz, y, alpha=0.6)
        axes_flat[ax_idx].set_xlabel('Feature 1')
        axes_flat[ax_idx].set_ylabel('Target')
    else:
        scatter = axes_flat[ax_idx].scatter(X_viz[:, 0], X_viz[:, 1], c=y, alpha=0.6, cmap='viridis')
        plt.colorbar(scatter, ax=axes_flat[ax_idx])
        axes_flat[ax_idx].set_xlabel('Feature 1')
        axes_flat[ax_idx].set_ylabel('Feature 2')
    
    axes_flat[ax_idx].set_title('Dados Originais')
    axes_flat[ax_idx].grid(True, alpha=0.3)
    
    # Plots dos modelos
    model_names = list(models.keys())
    for idx, model_name in enumerate(model_names[:4]):  # Mostrar até 4 modelos
        ax_idx = idx + 1  # +1 porque o primeiro subplot é para dados originais
        
        if ax_idx >= len(axes_flat):
            print(f"  Aviso: Não há espaço para mostrar {model_name}")
            continue
            
        if model_name in resultados_completos[cenario_nome]:
            resultados = resultados_completos[cenario_nome][model_name]
            
            if 'error' not in resultados:
                # CORREÇÃO: Criar scaler específico para visualização
                scaler_viz = StandardScaler()
                X_viz_scaled = scaler_viz.fit_transform(X_viz)
                
                # Criando superfície de predição para visualização
                if X_viz.shape[1] == 1:
                    x_plot = np.linspace(X_viz.min(), X_viz.max(), 100).reshape(-1, 1)
                    X_plot_scaled = scaler_viz.transform(x_plot)
                    y_plot = resultados['model'].predict(X_plot_scaled)
                    
                    axes_flat[ax_idx].scatter(X_viz, y, alpha=0.6, label='Dados')
                    axes_flat[ax_idx].plot(x_plot, y_plot, 'r-', linewidth=2, label='Predição')
                    axes_flat[ax_idx].set_xlabel('Feature 1')
                    
                else:
                    # Para 2D, criar meshgrid
                    x0 = np.linspace(X_viz[:, 0].min(), X_viz[:, 0].max(), 30)  # Reduzido para performance
                    x1 = np.linspace(X_viz[:, 1].min(), X_viz[:, 1].max(), 30)
                    X0, X1 = np.meshgrid(x0, x1)
                    X_plot = np.c_[X0.ravel(), X1.ravel()]
                    
                    # CORREÇÃO: Lidar corretamente com diferentes dimensionalidades
                    if X.shape[1] > 2:
                        # Para visualização 2D de dados multidimensionais
                        # Usar valores médios para as features extras
                        X_plot_full = np.zeros((X_plot.shape[0], X.shape[1]))
                        X_plot_full[:, :2] = X_plot
                        
                        # Preencher features extras com valores médios
                        for feat_idx in range(2, X.shape[1]):
                            X_plot_full[:, feat_idx] = np.mean(X[:, feat_idx])
                        
                        X_plot_scaled = scaler_cenario.transform(X_plot_full)
                    else:
                        X_plot_scaled = scaler_cenario.transform(X_plot)
                    
                    y_plot = resultados['model'].predict(X_plot_scaled)
                    y_plot = y_plot.reshape(X0.shape)
                    
                    # CORREÇÃO: Verificar se há dados válidos para contourf
                    if not np.all(np.isnan(y_plot)) and not np.all(np.isinf(y_plot)):
                        contour = axes_flat[ax_idx].contourf(X0, X1, y_plot, alpha=0.6, cmap='viridis', levels=20)
                        scatter_plot = axes_flat[ax_idx].scatter(X_viz[:, 0], X_viz[:, 1], c=y, alpha=0.8, cmap='viridis')
                        plt.colorbar(contour, ax=axes_flat[ax_idx])
                    else:
                        axes_flat[ax_idx].text(0.5, 0.5, 'Dados de predição\ninválidos', 
                                             ha='center', va='center', transform=axes_flat[ax_idx].transAxes)
                    
                    axes_flat[ax_idx].set_xlabel('Feature 1')
                    axes_flat[ax_idx].set_ylabel('Feature 2')
                
                axes_flat[ax_idx].set_title(f'{model_name}\nMSE: {resultados["mse_test"]:.4f}')
                axes_flat[ax_idx].grid(True, alpha=0.3)
                axes_flat[ax_idx].legend()
            else:
                # Mostrar mensagem de erro no subplot
                axes_flat[ax_idx].text(0.5, 0.5, f'Erro:\n{resultados["error"]}', 
                                     ha='center', va='center', transform=axes_flat[ax_idx].transAxes)
                axes_flat[ax_idx].set_title(f'{model_name} - Erro')
        else:
            # Modelo não foi treinado para este cenário
            axes_flat[ax_idx].text(0.5, 0.5, 'Modelo não\ntreinado', 
                                 ha='center', va='center', transform=axes_flat[ax_idx].transAxes)
            axes_flat[ax_idx].set_title(f'{model_name}')
    
    # Plot final: Comparação de performance (último subplot)
    ax_comp_idx = 5  # Último subplot em grid 2x3
    if ax_comp_idx < len(axes_flat):
        model_names_clean = []
        mse_scores = []
        
        for model_name in model_names:
            if (model_name in resultados_completos[cenario_nome] and 
                'error' not in resultados_completos[cenario_nome][model_name]):
                model_names_clean.append(model_name)
                mse_scores.append(resultados_completos[cenario_nome][model_name]['mse_test'])
        
        if model_names_clean:
            bars = axes_flat[ax_comp_idx].bar(model_names_clean, mse_scores, 
                                            color=plt.cm.Set3(np.arange(len(model_names_clean))))
            axes_flat[ax_comp_idx].set_title('Comparação de MSE (Teste)')
            axes_flat[ax_comp_idx].set_ylabel('MSE')
            axes_flat[ax_comp_idx].tick_params(axis='x', rotation=45)
            
            # Adicionar valores nas barras
            for bar, score in zip(bars, mse_scores):
                axes_flat[ax_comp_idx].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.001,
                            f'{score:.4f}', ha='center', va='bottom', fontsize=8)
        else:
            axes_flat[ax_comp_idx].text(0.5, 0.5, 'Sem dados\nválidos', 
                                      ha='center', va='center', transform=axes_flat[ax_comp_idx].transAxes)
            axes_flat[ax_comp_idx].set_title('Comparação de MSE (Teste)')
    
    # CORREÇÃO: Remover subplots não utilizados
    for i in range(len(axes_flat)):
        if i > ax_comp_idx and i < len(axes_flat):
            fig.delaxes(axes_flat[i])
    
    plt.tight_layout()
    plt.show()

'''
Análise de sensibilidade aos parâmetros - VERSÃO CORRIGIDA
'''
print("\n5. ANÁLISE DE SENSIBILIDADE")

# Focando no Kernel RBF para análise detalhada
print("\nAnálise de sensibilidade - Kernel RBF:")

for cenario_nome, X, y in cenarios:
    if 'KernelRidge RBF' in resultados_completos[cenario_nome]:
        resultados = resultados_completos[cenario_nome]['KernelRidge RBF']
        
        if 'error' not in resultados:
            print(f"\n{cenario_nome}:")
            print(f"  Melhores parâmetros: {resultados['best_params']}")
            print(f"  MSE treino: {resultados['mse_train']:.4f}")
            print(f"  MSE teste: {resultados['mse_test']:.4f}")
            print(f"  R² treino: {resultados['r2_train']:.4f}")
            print(f"  R² teste: {resultados['r2_test']:.4f}")

# CORREÇÃO: Adicionar análise de performance geral
print("\n6. ANÁLISE DE PERFORMANCE GERAL")

print("\nResumo dos Melhores Modelos por Cenário:")
for cenario_nome, resultados in resultados_completos.items():
    print(f"\n{cenario_nome}:")
    
    # Encontrar o melhor modelo baseado no MSE de teste
    melhor_modelo = None
    melhor_mse = float('inf')
    
    for model_name, dados in resultados.items():
        if 'error' not in dados and dados['mse_test'] < melhor_mse:
            melhor_mse = dados['mse_test']
            melhor_modelo = model_name
    
    if melhor_modelo:
        melhor_dados = resultados[melhor_modelo]
        print(f"  Melhor modelo: {melhor_modelo}")
        print(f"  MSE teste: {melhor_mse:.4f}")
        print(f"  R² teste: {melhor_dados['r2_test']:.4f}")
        print(f"  Parâmetros: {melhor_dados['best_params']}")

'''
Recomendações práticas - VERSÃO ATUALIZADA
'''
print("\n7. RECOMENDAÇÕES PRÁTICAS")

print("\nBaseado na análise experimental:")

print("\nEscolha do Kernel:")
print("  - Linear: Para relações aproximadamente lineares (MSE: 0.0483 no cenário linear)")
print("  - RBF: Para relações não-lineares suaves (MSE: 0.0359 no Friedman1)")
print("  - Polinomial: Para relações polinomiais conhecidas (MSE: 0.0073 no Friedman1)")

print("\nOtimização de Hiperparâmetros:")
print("  - Sempre usar validação cruzada para alpha e gamma")
print("  - Gamma controla a complexidade do modelo (valores típicos: 0.1-10)")
print("  - Alpha controla a regularização (valores típicos: 0.1-100)")

print("\nCenários de Aplicação:")
print("  - KernelRidge é preferível quando há relações não-lineares")
print("  - Ridge linear é suficiente para dados linearmente separáveis")
print("  - Kernel polinomial performou excepcionalmente bem no Friedman1")

print("\nVantagens do KernelRidge:")
print("  + Flexibilidade para capturar não-linearidades")
print("  + Solução em forma fechada")
print("  + Boa generalização com tuning adequado")
print("  - Custo computacional mais alto")
print("  - Menos interpretável que modelos lineares")

print("\nObservações dos Resultados:")
print("  - No cenário complexo, KernelRidge polinomial teve R² de 0.9045 vs 0.0287 do Ridge linear")
print("  - KernelRidge linear performou similarmente ao Ridge linear tradicional")
print("  - O tuning de parâmetros é crucial para o desempenho")

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

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

import numpy as np

import matplotlib.pyplot as plt

from sklearn.kernel_ridge import KernelRidge

from sklearn.linear_model import Ridge

from sklearn.model_selection import train_test_split, GridSearchCV

from sklearn.metrics import mean_squared_error, r2_score

from sklearn.preprocessing import StandardScaler

from sklearn.datasets import make_regression, make_friedman1

import warnings

warnings.filterwarnings('ignore')

'''

Estudo comparativo da Regressão de Crista do Kernel

em diferentes cenários de dados

'''

print("=== REGRESSÃO DE CRISTA DO KERNEL ===")

# Criando diferentes cenários de dados

print("\n1. CONFIGURAÇÃO DOS CENÁRIOS")

# Cenário 1: Relação linear com ruído

X_linear, y_linear = make_regression(

n_samples=200, n_features=2, noise=10, random_state=42

)

# Cenário 2: Relação não-linear (Friedman1)

X_nonlinear, y_nonlinear = make_friedman1(

n_samples=200, n_features=5, noise=0.1, random_state=42

)

# Cenário 3: Dados com interações complexas

np.random.seed(42)

X_complex = np.random.randn(200, 3)

y_complex = (X_complex[:, 0]**2 +

np.sin(X_complex[:, 1] * 2) +

X_complex[:, 0] * X_complex[:, 2] +

np.random.normal(0, 0.1, 200))

cenarios = [

('Linear', X_linear, y_linear),

('Não-Linear (Friedman1)', X_nonlinear, y_nonlinear),

('Complexo', X_complex, y_complex)

]

'''

Configuração dos modelos comparativos

'''

print("\n2. CONFIGURAÇÃO DOS MODELOS")

models = {

'Ridge Linear': Ridge(),

'KernelRidge RBF': KernelRidge(kernel='rbf'),

'KernelRidge Polinomial': KernelRidge(kernel='polynomial'),

'KernelRidge Linear': KernelRidge(kernel='linear')

}

# Espaço de parâmetros para tuning

param_grids = {

'Ridge Linear': {'alpha': [0.1, 1.0, 10.0, 100.0]},

'KernelRidge RBF': {

'alpha': [0.1, 1.0, 10.0],

'gamma': [0.1, 1.0, 10.0]

'KernelRidge Polinomial': {

'alpha': [0.1, 1.0, 10.0],

'degree': [2, 3, 4],

'gamma': [0.1, 1.0]

'KernelRidge Linear': {

'alpha': [0.1, 1.0, 10.0, 100.0]

}

'''

Avaliação sistemática em todos os cenários

'''

print("\n3. AVALIAÇÃO COMPARATIVA")

resultados_completos = {}

for cenario_nome, X, y in cenarios:

print(f"\n--- {cenario_nome} ---")

print(f"Dimensões: {X.shape}")

# Padronizando os dados

scaler = StandardScaler()

X_scaled = scaler.fit_transform(X)

y_scaled = StandardScaler().fit_transform(y.reshape(-1, 1)).ravel()

# Dividindo em treino e teste

X_train, X_test, y_train, y_test = train_test_split(

X_scaled, y_scaled, test_size=0.3, random_state=42

)

cenario_resultados = {}

for model_name in models.keys():

print(f" Treinando {model_name}...")

try:

# Grid Search com validação cruzada

grid_search = GridSearchCV(

models[model_name],

param_grids[model_name],

cv=5,

scoring='neg_mean_squared_error',

n_jobs=-1

)

grid_search.fit(X_train, y_train)

# Melhor modelo

best_model = grid_search.best_estimator_

best_params = grid_search.best_params_

# Previsões e métricas

y_pred_train = best_model.predict(X_train)

y_pred_test = best_model.predict(X_test)

mse_train = mean_squared_error(y_train, y_pred_train)

mse_test = mean_squared_error(y_test, y_pred_test)

r2_train = r2_score(y_train, y_pred_train)

r2_test = r2_score(y_test, y_pred_test)

# Armazenando resultados

cenario_resultados[model_name] = {

'model': best_model,

'best_params': best_params,

'mse_train': mse_train,

'mse_test': mse_test,

'r2_train': r2_train,

'r2_test': r2_test,

'best_score': -grid_search.best_score_ # Convertendo para MSE positivo

}

print(f" Melhores parâmetros: {best_params}")

print(f" MSE teste: {mse_test:.4f}, R² teste: {r2_test:.4f}")

except Exception as e:

print(f" ERRO: {e}")

cenario_resultados[model_name] = {'error': str(e)}

resultados_completos[cenario_nome] = cenario_resultados

'''

Visualização comparativa dos resultados - VERSÃO CORRIGIDA

'''

print("\n4. VISUALIZAÇÃO DOS RESULTADOS")

# Criando visualizações para cada cenário

for cenario_idx, (cenario_nome, X, y) in enumerate(cenarios):

# CORREÇÃO: Usar o scaler específico para cada cenário

scaler_cenario = StandardScaler()

X_scaled_cenario = scaler_cenario.fit_transform(X)

y_scaled_cenario = StandardScaler().fit_transform(y.reshape(-1, 1)).ravel()

# Para visualização, usar apenas os dois primeiros features se houver mais

if X.shape[1] > 2:

X_viz = X[:, :2]

print(f" {cenario_nome}: Usando 2 features para visualização de {X.shape[1]} features totais")

else:

X_viz = X

# CORREÇÃO: Ajustar layout dos subplots baseado no número de modelos

n_models = min(4, len(models)) # Máximo 4 modelos para visualização

n_rows = 2

n_cols = 3 # 1 para dados originais + 4 para modelos + 1 para comparação

fig, axes = plt.subplots(n_rows, n_cols, figsize=(18, 10))

fig.suptitle(f'Regressão de Crista do Kernel - {cenario_nome}', fontsize=16)

# Aplanar axes para facilitar acesso

if n_rows > 1 and n_cols > 1:

axes_flat = axes.flatten()

else:

axes_flat = [axes] if n_cols == 1 else axes

# Plot 1: Dados originais

ax_idx = 0

if X_viz.shape[1] == 1:

axes_flat[ax_idx].scatter(X_viz, y, alpha=0.6)

axes_flat[ax_idx].set_xlabel('Feature 1')

axes_flat[ax_idx].set_ylabel('Target')

else:

scatter = axes_flat[ax_idx].scatter(X_viz[:, 0], X_viz[:, 1], c=y, alpha=0.6, cmap='viridis')

plt.colorbar(scatter, ax=axes_flat[ax_idx])

axes_flat[ax_idx].set_xlabel('Feature 1')

axes_flat[ax_idx].set_ylabel('Feature 2')

axes_flat[ax_idx].set_title('Dados Originais')

axes_flat[ax_idx].grid(True, alpha=0.3)

# Plots dos modelos

model_names = list(models.keys())

for idx, model_name in enumerate(model_names[:4]): # Mostrar até 4 modelos

ax_idx = idx + 1 # +1 porque o primeiro subplot é para dados originais

if ax_idx >= len(axes_flat):

print(f" Aviso: Não há espaço para mostrar {model_name}")

continue

if model_name in resultados_completos[cenario_nome]:

resultados = resultados_completos[cenario_nome][model_name]

if 'error' not in resultados:

# CORREÇÃO: Criar scaler específico para visualização

scaler_viz = StandardScaler()

X_viz_scaled = scaler_viz.fit_transform(X_viz)

# Criando superfície de predição para visualização

if X_viz.shape[1] == 1:

x_plot = np.linspace(X_viz.min(), X_viz.max(), 100).reshape(-1, 1)

X_plot_scaled = scaler_viz.transform(x_plot)

y_plot = resultados['model'].predict(X_plot_scaled)

axes_flat[ax_idx].scatter(X_viz, y, alpha=0.6, label='Dados')

axes_flat[ax_idx].plot(x_plot, y_plot, 'r-', linewidth=2, label='Predição')

axes_flat[ax_idx].set_xlabel('Feature 1')

else:

# Para 2D, criar meshgrid

x0 = np.linspace(X_viz[:, 0].min(), X_viz[:, 0].max(), 30) # Reduzido para performance

x1 = np.linspace(X_viz[:, 1].min(), X_viz[:, 1].max(), 30)

X0, X1 = np.meshgrid(x0, x1)

X_plot = np.c_[X0.ravel(), X1.ravel()]

# CORREÇÃO: Lidar corretamente com diferentes dimensionalidades

if X.shape[1] > 2:

# Para visualização 2D de dados multidimensionais

# Usar valores médios para as features extras

X_plot_full = np.zeros((X_plot.shape[0], X.shape[1]))

X_plot_full[:, :2] = X_plot

# Preencher features extras com valores médios

for feat_idx in range(2, X.shape[1]):

X_plot_full[:, feat_idx] = np.mean(X[:, feat_idx])

X_plot_scaled = scaler_cenario.transform(X_plot_full)

else:

X_plot_scaled = scaler_cenario.transform(X_plot)

y_plot = resultados['model'].predict(X_plot_scaled)

y_plot = y_plot.reshape(X0.shape)

# CORREÇÃO: Verificar se há dados válidos para contourf

if not np.all(np.isnan(y_plot)) and not np.all(np.isinf(y_plot)):

contour = axes_flat[ax_idx].contourf(X0, X1, y_plot, alpha=0.6, cmap='viridis', levels=20)

scatter_plot = axes_flat[ax_idx].scatter(X_viz[:, 0], X_viz[:, 1], c=y, alpha=0.8, cmap='viridis')

plt.colorbar(contour, ax=axes_flat[ax_idx])

else:

axes_flat[ax_idx].text(0.5, 0.5, 'Dados de predição\ninválidos',

ha='center', va='center', transform=axes_flat[ax_idx].transAxes)

axes_flat[ax_idx].set_xlabel('Feature 1')

axes_flat[ax_idx].set_ylabel('Feature 2')

axes_flat[ax_idx].set_title(f'{model_name}\nMSE: {resultados["mse_test"]:.4f}')

axes_flat[ax_idx].grid(True, alpha=0.3)

axes_flat[ax_idx].legend()

else:

# Mostrar mensagem de erro no subplot

axes_flat[ax_idx].text(0.5, 0.5, f'Erro:\n{resultados["error"]}',

ha='center', va='center', transform=axes_flat[ax_idx].transAxes)

axes_flat[ax_idx].set_title(f'{model_name} - Erro')

else:

# Modelo não foi treinado para este cenário

axes_flat[ax_idx].text(0.5, 0.5, 'Modelo não\ntreinado',

ha='center', va='center', transform=axes_flat[ax_idx].transAxes)

axes_flat[ax_idx].set_title(f'{model_name}')

# Plot final: Comparação de performance (último subplot)

ax_comp_idx = 5 # Último subplot em grid 2x3

if ax_comp_idx < len(axes_flat):

model_names_clean = []

mse_scores = []

for model_name in model_names:

if (model_name in resultados_completos[cenario_nome] and

'error' not in resultados_completos[cenario_nome][model_name]):

model_names_clean.append(model_name)

mse_scores.append(resultados_completos[cenario_nome][model_name]['mse_test'])

if model_names_clean:

bars = axes_flat[ax_comp_idx].bar(model_names_clean, mse_scores,

color=plt.cm.Set3(np.arange(len(model_names_clean))))

axes_flat[ax_comp_idx].set_title('Comparação de MSE (Teste)')

axes_flat[ax_comp_idx].set_ylabel('MSE')

axes_flat[ax_comp_idx].tick_params(axis='x', rotation=45)

# Adicionar valores nas barras

for bar, score in zip(bars, mse_scores):

axes_flat[ax_comp_idx].text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.001,

f'{score:.4f}', ha='center', va='bottom', fontsize=8)

else:

axes_flat[ax_comp_idx].text(0.5, 0.5, 'Sem dados\nválidos',

ha='center', va='center', transform=axes_flat[ax_comp_idx].transAxes)

axes_flat[ax_comp_idx].set_title('Comparação de MSE (Teste)')

# CORREÇÃO: Remover subplots não utilizados

for i in range(len(axes_flat)):

if i > ax_comp_idx and i < len(axes_flat):

fig.delaxes(axes_flat[i])

plt.tight_layout()

plt.show()

'''

Análise de sensibilidade aos parâmetros - VERSÃO CORRIGIDA

'''

print("\n5. ANÁLISE DE SENSIBILIDADE")

# Focando no Kernel RBF para análise detalhada

print("\nAnálise de sensibilidade - Kernel RBF:")

for cenario_nome, X, y in cenarios:

if 'KernelRidge RBF' in resultados_completos[cenario_nome]:

resultados = resultados_completos[cenario_nome]['KernelRidge RBF']

if 'error' not in resultados:

print(f"\n{cenario_nome}:")

print(f" Melhores parâmetros: {resultados['best_params']}")

print(f" MSE treino: {resultados['mse_train']:.4f}")

print(f" MSE teste: {resultados['mse_test']:.4f}")

print(f" R² treino: {resultados['r2_train']:.4f}")

print(f" R² teste: {resultados['r2_test']:.4f}")

# CORREÇÃO: Adicionar análise de performance geral

print("\n6. ANÁLISE DE PERFORMANCE GERAL")

print("\nResumo dos Melhores Modelos por Cenário:")

for cenario_nome, resultados in resultados_completos.items():

print(f"\n{cenario_nome}:")

# Encontrar o melhor modelo baseado no MSE de teste

melhor_modelo = None

melhor_mse = float('inf')

for model_name, dados in resultados.items():

if 'error' not in dados and dados['mse_test'] < melhor_mse:

melhor_mse = dados['mse_test']

melhor_modelo = model_name

if melhor_modelo:

melhor_dados = resultados[melhor_modelo]

print(f" Melhor modelo: {melhor_modelo}")

print(f" MSE teste: {melhor_mse:.4f}")

print(f" R² teste: {melhor_dados['r2_test']:.4f}")

print(f" Parâmetros: {melhor_dados['best_params']}")

'''

Recomendações práticas - VERSÃO ATUALIZADA

'''

print("\n7. RECOMENDAÇÕES PRÁTICAS")

print("\nBaseado na análise experimental:")

print("\nEscolha do Kernel:")

print(" - Linear: Para relações aproximadamente lineares (MSE: 0.0483 no cenário linear)")

print(" - RBF: Para relações não-lineares suaves (MSE: 0.0359 no Friedman1)")

print(" - Polinomial: Para relações polinomiais conhecidas (MSE: 0.0073 no Friedman1)")

print("\nOtimização de Hiperparâmetros:")

print(" - Sempre usar validação cruzada para alpha e gamma")

print(" - Gamma controla a complexidade do modelo (valores típicos: 0.1-10)")

print(" - Alpha controla a regularização (valores típicos: 0.1-100)")

print("\nCenários de Aplicação:")

print(" - KernelRidge é preferível quando há relações não-lineares")

print(" - Ridge linear é suficiente para dados linearmente separáveis")

print(" - Kernel polinomial performou excepcionalmente bem no Friedman1")

print("\nVantagens do KernelRidge:")

print(" + Flexibilidade para capturar não-linearidades")

print(" + Solução em forma fechada")

print(" + Boa generalização com tuning adequado")

print(" - Custo computacional mais alto")

print(" - Menos interpretável que modelos lineares")

print("\nObservações dos Resultados:")

print(" - No cenário complexo, KernelRidge polinomial teve R² de 0.9045 vs 0.0287 do Ridge linear")

print(" - KernelRidge linear performou similarmente ao Ridge linear tradicional")

print(" - O tuning de parâmetros é crucial para o desempenho")

Interpretação dos Resultados

Analisando os experimentos comparativos, podemos observar padrões importantes:

O KernelRidge com kernel RBF geralmente performa melhor em dados não-lineares
A regressão Ridge linear mantém vantagem em dados verdadeiramente lineares
O tuning de hiperparâmetros é crucial para o desempenho do KernelRidge
O custo computacional aumenta significativamente com kernels não-lineares

Considerações de Implementação

Complexidade Computacional

A Regressão de Crista do Kernel possui complexidade O(n³) para o treinamento devido à inversão da matriz do kernel. Ademais, o armazenamento da matriz do kernel requer O(n²) de memória.

Escalabilidade

Para datasets grandes, estratégias alternativas podem ser necessárias:

Uso de aproximações de kernel (Nyström method)
Amostragem para reduzir o tamanho do dataset
Implementações distribuídas

Conexões com o Ecossistema scikit-learn

Similarmente a outros estimadores do scikit-learn, o KernelRidge oferece:

Interface consistente com fit, predict, e score
Suporte a GridSearchCV para tuning de hiperparâmetros
Integração com pipelines de pré-processamento
Compatibilidade com métricas de avaliação padrão

Conclusão

A Regressão de Crista do Kernel representa uma extensão poderosa dos métodos de regularização linear, incorporando a flexibilidade dos kernels para lidar com relações complexas nos dados. Embora computacionalmente mais custosa que abordagens lineares, oferece capacidade de modelagem significativamente expandida.

Portanto, esta técnica é particularmente valiosa em cenários onde relações não-lineares estão presentes e onde o custo computacional é aceitável, preenchendo uma importante lacuna entre modelos lineares simples e métodos não-paramétricos complexos.

Referência

Este post explora o item 1.3. Regressão de Crista do Kernel da documentação do scikit-learn:

https://scikit-learn.org/0.21/modules/kernel_ridge.html

Lasso LARS e a fusão entre Regularização L1 e Abordagem Geométric

19/12/202517/10/2025 Por antonino

Introdução ao Lasso LARS

O Lasso LARS constitui uma implementação eficiente que combina as vantagens da regressão Lasso com o algoritmo Least Angle Regression. Primordialmente, esta abordagem permite computar a solução completa do caminho Lasso de forma computacionalmente eficiente, sem a necessidade de otimização convexa iterativa.

Princípio de Funcionamento

O Lasso LARS opera através de um processo sequencial similar ao LARS tradicional, mas com uma modificação crucial: quando um coeficiente atinge zero, a feature correspondente é removida do conjunto ativo. Esta simples alteração garante que as soluções produzidas sejam idênticas às do Lasso convencional.

Diferenciação do LARS Clássico

Enquanto o LARS tradicional sempre mantém as features no modelo uma vez selecionadas, o Lasso LARS permite a remoção de features quando seus coeficientes se tornam zero. Isto produz modelos mais esparsos e alinhados com a filosofia de seleção de features do Lasso.

Vantagens Computacionais

Computação eficiente de todo o caminho de regularização
Menor custo computacional para problemas com muitas features
Estabilidade numérica superior em comparação com métodos baseados em descida de gradiente
Capacidade de lidar com situações onde número de features excede número de amostras

Formulação Matemática

O Lasso LARS resolve o mesmo problema de otimização do Lasso tradicional:

\(\min_{w} \frac{1}{2n_{\text{samples}}} ||X w – y||_2^2 + \alpha ||w||_1\)

Contudo, a estratégia de solução é fundamentalmente diferente. Ao invés de otimização convexa, utiliza-se uma abordagem geométrica baseada em direções equiangulares.

Algoritmo Lasso LARS

Iniciar com coeficientes zero e resíduo igual ao target
Encontrar a feature mais correlacionada com o resíduo
Mover na direção desta feature até que outra feature tenha correlação igual
Continuar na direção equiangular entre features ativas
Se algum coeficiente atingir zero, remover a feature do conjunto ativo
Repetir até inclusão de todas as features relevantes

Cenários de Aplicação Ideal

O Lasso LARS é particularmente adequado para:

Problemas com número de features maior que número de amostras (p > n)
Análise exploratória onde se deseja examinar todo o caminho de regularização
Situações que requerem seleção eficiente de features
Aplicações onde a estabilidade numérica é crítica

Implementação no scikit-learn

No scikit-learn, a classe LassoLars implementa esta técnica. Ademais, LassoLarsCV fornece seleção automática do parâmetro alpha através de validação cruzada.

Parâmetros Principais

alpha: Parâmetro de regularização (constante ou array)
fit_intercept: Cálculo do termo intercepto
verbose: Nível de detalhe durante execução
normalize: Normalização prévia dos dados
max_iter: Número máximo de iterações

Exemplo Prático: Comparação com Outras Abordagens

O exemplo a seguir demonstra o uso do Lasso LARS e sua comparação com implementações alternativas do Lasso:

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LassoLars, Lasso, LassoLarsIC
from sklearn.datasets import make_regression
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error, r2_score
from sklearn.preprocessing import StandardScaler
import time

print("=" * 60)
print("COMPARAÇÃO: LASSO LARS VS LASSO TRADICIONAL")
print("=" * 60)

# Gerar dados com alta dimensionalidade (p > n scenario)
X, y = make_regression(n_samples=80, n_features=50, n_informative=10, 
                       noise=0.3, random_state=42)

# Normalizar os dados
scaler = StandardScaler()
X_scaled = scaler.fit_transform(X)
y_centered = y - np.mean(y)

# Dividir em treino e teste
X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_centered, 
                                                    test_size=0.3, random_state=42)

print(f"Cenário p > n: {X.shape[1]} features, {X.shape[0]} amostras")
print(f"Features informativas verdadeiras: 10")

# 1. Lasso LARS com alpha fixo
print("\n1. LASSO LARS (ALPHA=0.1)")
start_time = time.time()
lasso_lars = LassoLars(alpha=0.1)
lasso_lars.fit(X_train, y_train)
lars_time = time.time() - start_time

y_pred_lars = lasso_lars.predict(X_test)
coef_nao_nulos_lars = np.sum(lasso_lars.coef_ != 0)

print(f"Tempo de execução: {lars_time:.4f} segundos")
print(f"Coeficientes não nulos: {coef_nao_nulos_lars}/50")
print(f"MSE: {mean_squared_error(y_test, y_pred_lars):.4f}")
print(f"R²: {r2_score(y_test, y_pred_lars):.4f}")

# 2. Lasso tradicional
print("\n2. LASSO TRADICIONAL (ALPHA=0.1)")
start_time = time.time()
lasso_trad = Lasso(alpha=0.1, max_iter=10000, random_state=42)
lasso_trad.fit(X_train, y_train)
lasso_time = time.time() - start_time

y_pred_trad = lasso_trad.predict(X_test)
coef_nao_nulos_trad = np.sum(lasso_trad.coef_ != 0)

print(f"Tempo de execução: {lasso_time:.4f} segundos")
print(f"Coeficientes não nulos: {coef_nao_nulos_trad}/50")
print(f"MSE: {mean_squared_error(y_test, y_pred_trad):.4f}")
print(f"R²: {r2_score(y_test, y_pred_trad):.4f}")

# 3. Lasso LARS com critério de informação
print("\n3. LASSO LARS COM CRITÉRIO AIC")
lasso_lars_aic = LassoLarsIC(criterion='aic')
lasso_lars_aic.fit(X_train, y_train)
y_pred_aic = lasso_lars_aic.predict(X_test)

print(f"Alpha selecionado (AIC): {lasso_lars_aic.alpha_:.6f}")
print(f"Coeficientes não nulos: {np.sum(lasso_lars_aic.coef_ != 0)}/50")
print(f"MSE: {mean_squared_error(y_test, y_pred_aic):.4f}")

# Visualização comparativa
plt.figure(figsize=(15, 10))

# Gráfico 1: Comparação de coeficientes
plt.subplot(2, 2, 1)
coeficientes = [lasso_lars.coef_, lasso_trad.coef_, lasso_lars_aic.coef_]
nomes = ['LassoLars', 'Lasso Trad', 'LassoLars AIC']
cores = ['blue', 'red', 'green']

for i, (coef, nome) in enumerate(zip(coeficientes, nomes)):
    plt.bar(np.arange(50) + i*0.25, coef, width=0.23, 
            label=nome, alpha=0.7, color=cores[i])

plt.axhline(y=0, color='black', linestyle='-', alpha=0.3)
plt.xlabel('Features')
plt.ylabel('Valor do Coeficiente')
plt.title('Comparação dos Coeficientes - Cenário p > n')
plt.legend()
plt.grid(True, alpha=0.3)

# Gráfico 2: Caminho de regularização LassoLars
plt.subplot(2, 2, 2)
# Calcular caminho completo com alpha zero
lars_path = LassoLars(alpha=0)
lars_path.fit(X_train, y_train)

# Extrair alphas ao longo do caminho
alphas = lars_path.alphas_
coef_path = lars_path.coef_path_

# Plotar caminho para primeiras 15 features
for i in range(min(15, coef_path.shape[0])):
    plt.plot(alphas, coef_path[i], label=f'F{i+1}', linewidth=1.5)

plt.xlabel('Alpha')
plt.ylabel('Valor do Coeficiente')
plt.title('Caminho LassoLars - Evolução com Alpha')
plt.xscale('log')
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')
plt.grid(True, alpha=0.3)

# Gráfico 3: Comparação de esparsidade
plt.subplot(2, 2, 3)
features_nao_nulas = [np.sum(coef != 0) for coef in coeficientes]
plt.bar(nomes, features_nao_nulas, color=cores, alpha=0.7)
plt.ylabel('Número de Features Selecionadas')
plt.title('Comparação de Esparsidade')
for i, v in enumerate(features_nao_nulas):
    plt.text(i, v + 0.5, str(v), ha='center', va='bottom')
plt.grid(True, alpha=0.3)

# Gráfico 4: Comparação de performance
plt.subplot(2, 2, 4)
mse_valores = [
    mean_squared_error(y_test, y_pred_lars),
    mean_squared_error(y_test, y_pred_trad),
    mean_squared_error(y_test, y_pred_aic)
]
tempos = [lars_time, lasso_time, 0]  # AIC não tem tempo separado

x_pos = np.arange(len(nomes))
width = 0.35

plt.bar(x_pos - width/2, mse_valores, width, label='MSE', alpha=0.7, color='blue')
plt.bar(x_pos + width/2, tempos, width, label='Tempo (s)', alpha=0.7, color='orange')
plt.xlabel('Método')
plt.ylabel('MSE / Tempo (s)')
plt.title('Performance e Eficiência Computacional')
plt.xticks(x_pos, nomes)
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Análise de concordância na seleção de features
print("\n4. ANÁLISE DE CONCORDÂNCIA NA SELEÇÃO")
print("Feature | LassoLars | LassoTrad | Concordância")
print("-" * 50)
concordancia_total = 0
for i in range(50):
    selecionada_lars = lasso_lars.coef_[i] != 0
    selecionada_trad = lasso_trad.coef_[i] != 0
    concordancia = "SIM" if selecionada_lars == selecionada_trad else "NÃO"
    if selecionada_lars == selecionada_trad:
        concordancia_total += 1
    
    status_lars = "SEL" if selecionada_lars else "---"
    status_trad = "SEL" if selecionada_trad else "---"
    print(f"F{i+1:2d}     |    {status_lars}    |    {status_trad}    |     {concordancia}")

print(f"\nTaxa de concordância: {concordancia_total/50*100:.1f}%")

# Análise adicional do caminho LassoLars
print("\n5. INFORMAÇÕES DO CAMINHO LASSO LARS")
print(f"Número de iterações: {len(lars_path.alphas_)}")
print(f"Alpha máximo: {lars_path.alphas_[0]:.6f}")
print(f"Alpha mínimo: {lars_path.alphas_[-1]:.6f}")

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

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

import numpy as np

import matplotlib.pyplot as plt

from sklearn.linear_model import LassoLars, Lasso, LassoLarsIC

from sklearn.datasets import make_regression

from sklearn.model_selection import train_test_split

from sklearn.metrics import mean_squared_error, r2_score

from sklearn.preprocessing import StandardScaler

import time

print("=" * 60)

print("COMPARAÇÃO: LASSO LARS VS LASSO TRADICIONAL")

print("=" * 60)

# Gerar dados com alta dimensionalidade (p > n scenario)

X, y = make_regression(n_samples=80, n_features=50, n_informative=10,

noise=0.3, random_state=42)

# Normalizar os dados

scaler = StandardScaler()

X_scaled = scaler.fit_transform(X)

y_centered = y - np.mean(y)

# Dividir em treino e teste

X_train, X_test, y_train, y_test = train_test_split(X_scaled, y_centered,

test_size=0.3, random_state=42)

print(f"Cenário p > n: {X.shape[1]} features, {X.shape[0]} amostras")

print(f"Features informativas verdadeiras: 10")

# 1. Lasso LARS com alpha fixo

print("\n1. LASSO LARS (ALPHA=0.1)")

start_time = time.time()

lasso_lars = LassoLars(alpha=0.1)

lasso_lars.fit(X_train, y_train)

lars_time = time.time() - start_time

y_pred_lars = lasso_lars.predict(X_test)

coef_nao_nulos_lars = np.sum(lasso_lars.coef_ != 0)

print(f"Tempo de execução: {lars_time:.4f} segundos")

print(f"Coeficientes não nulos: {coef_nao_nulos_lars}/50")

print(f"MSE: {mean_squared_error(y_test, y_pred_lars):.4f}")

print(f"R²: {r2_score(y_test, y_pred_lars):.4f}")

# 2. Lasso tradicional

print("\n2. LASSO TRADICIONAL (ALPHA=0.1)")

start_time = time.time()

lasso_trad = Lasso(alpha=0.1, max_iter=10000, random_state=42)

lasso_trad.fit(X_train, y_train)

lasso_time = time.time() - start_time

y_pred_trad = lasso_trad.predict(X_test)

coef_nao_nulos_trad = np.sum(lasso_trad.coef_ != 0)

print(f"Tempo de execução: {lasso_time:.4f} segundos")

print(f"Coeficientes não nulos: {coef_nao_nulos_trad}/50")

print(f"MSE: {mean_squared_error(y_test, y_pred_trad):.4f}")

print(f"R²: {r2_score(y_test, y_pred_trad):.4f}")

# 3. Lasso LARS com critério de informação

print("\n3. LASSO LARS COM CRITÉRIO AIC")

lasso_lars_aic = LassoLarsIC(criterion='aic')

lasso_lars_aic.fit(X_train, y_train)

y_pred_aic = lasso_lars_aic.predict(X_test)

print(f"Alpha selecionado (AIC): {lasso_lars_aic.alpha_:.6f}")

print(f"Coeficientes não nulos: {np.sum(lasso_lars_aic.coef_ != 0)}/50")

print(f"MSE: {mean_squared_error(y_test, y_pred_aic):.4f}")

# Visualização comparativa

plt.figure(figsize=(15, 10))

# Gráfico 1: Comparação de coeficientes

plt.subplot(2, 2, 1)

coeficientes = [lasso_lars.coef_, lasso_trad.coef_, lasso_lars_aic.coef_]

nomes = ['LassoLars', 'Lasso Trad', 'LassoLars AIC']

cores = ['blue', 'red', 'green']

for i, (coef, nome) in enumerate(zip(coeficientes, nomes)):

plt.bar(np.arange(50) + i*0.25, coef, width=0.23,

label=nome, alpha=0.7, color=cores[i])

plt.axhline(y=0, color='black', linestyle='-', alpha=0.3)

plt.xlabel('Features')

plt.ylabel('Valor do Coeficiente')

plt.title('Comparação dos Coeficientes - Cenário p > n')

plt.legend()

plt.grid(True, alpha=0.3)

# Gráfico 2: Caminho de regularização LassoLars

plt.subplot(2, 2, 2)

# Calcular caminho completo com alpha zero

lars_path = LassoLars(alpha=0)

lars_path.fit(X_train, y_train)

# Extrair alphas ao longo do caminho

alphas = lars_path.alphas_

coef_path = lars_path.coef_path_

# Plotar caminho para primeiras 15 features

for i in range(min(15, coef_path.shape[0])):

plt.plot(alphas, coef_path[i], label=f'F{i+1}', linewidth=1.5)

plt.xlabel('Alpha')

plt.ylabel('Valor do Coeficiente')

plt.title('Caminho LassoLars - Evolução com Alpha')

plt.xscale('log')

plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left')

plt.grid(True, alpha=0.3)

# Gráfico 3: Comparação de esparsidade

plt.subplot(2, 2, 3)

features_nao_nulas = [np.sum(coef != 0) for coef in coeficientes]

plt.bar(nomes, features_nao_nulas, color=cores, alpha=0.7)

plt.ylabel('Número de Features Selecionadas')

plt.title('Comparação de Esparsidade')

for i, v in enumerate(features_nao_nulas):

plt.text(i, v + 0.5, str(v), ha='center', va='bottom')

plt.grid(True, alpha=0.3)

# Gráfico 4: Comparação de performance

plt.subplot(2, 2, 4)

mse_valores = [

mean_squared_error(y_test, y_pred_lars),

mean_squared_error(y_test, y_pred_trad),

mean_squared_error(y_test, y_pred_aic)

]

tempos = [lars_time, lasso_time, 0] # AIC não tem tempo separado

x_pos = np.arange(len(nomes))

width = 0.35

plt.bar(x_pos - width/2, mse_valores, width, label='MSE', alpha=0.7, color='blue')

plt.bar(x_pos + width/2, tempos, width, label='Tempo (s)', alpha=0.7, color='orange')

plt.xlabel('Método')

plt.ylabel('MSE / Tempo (s)')

plt.title('Performance e Eficiência Computacional')

plt.xticks(x_pos, nomes)

plt.legend()

plt.grid(True, alpha=0.3)

plt.tight_layout()

plt.show()

# Análise de concordância na seleção de features

print("\n4. ANÁLISE DE CONCORDÂNCIA NA SELEÇÃO")

print("Feature | LassoLars | LassoTrad | Concordância")

print("-" * 50)

concordancia_total = 0

for i in range(50):

selecionada_lars = lasso_lars.coef_[i] != 0

selecionada_trad = lasso_trad.coef_[i] != 0

concordancia = "SIM" if selecionada_lars == selecionada_trad else "NÃO"

if selecionada_lars == selecionada_trad:

concordancia_total += 1

status_lars = "SEL" if selecionada_lars else "---"

status_trad = "SEL" if selecionada_trad else "---"

print(f"F{i+1:2d} | {status_lars} | {status_trad} | {concordancia}")

print(f"\nTaxa de concordância: {concordancia_total/50*100:.1f}%")

# Análise adicional do caminho LassoLars

print("\n5. INFORMAÇÕES DO CAMINHO LASSO LARS")

print(f"Número de iterações: {len(lars_path.alphas_)}")

print(f"Alpha máximo: {lars_path.alphas_[0]:.6f}")

print(f"Alpha mínimo: {lars_path.alphas_[-1]:.6f}")

Considerações sobre Estabilidade Numérica

O Lasso LARS demonstra particular robustez em situações numericamente desafiadoras. Por utilizar decomposições matriciais estáveis e evitar problemas de convergência associados a métodos iterativos, frequentemente produz resultados mais confiáveis em problemas mal condicionados.

Limitações e Considerações

Embora poderoso, o Lasso LARS apresenta algumas limitações:

Pode ser menos eficiente quando número de amostras é muito grande
Implementação atual não suporta sparse matrices de forma nativa
Menos flexível para personalização em comparação com abordagens baseadas em otimização

Considerações Finais

Inegavelmente, o Lasso LARS representa uma fusão elegante entre a intuitividade geométrica do LARS e o poder de seleção de features do Lasso. Sua capacidade de computar eficientemente todo o caminho de regularização o torna invaluable para análise exploratória e seleção de modelos.

Decerto, em cenários de alta dimensionalidade (p > n), o Lasso LARS frequentemente supera implementações tradicionais do Lasso tanto em termos de eficiência computacional quanto de estabilidade numérica. Ademais, fornece insights valiosos sobre a ordem de importância das features através do exame do caminho de regularização.

Analogamente a outras técnicas de regularização, a escolha do parâmetro alpha permanece crucial. Portanto, a combinação do Lasso LARS com métodos de validação cruzada ou critérios de informação constitui uma abordagem robusta para problemas práticos de machine learning.